UNCLASSIFIED 


_ AD  NUMBER _ 

AD441983 

LIMITATION  CHANGES 

TO: 

Approved  for  public  release;  distribution  is 
unlimited. 


FROM: 

Distribution  authorized  to  U.S.  Gov't,  agencies 
and  their  contractors; 

Administrative/Operational  Use;  23  APR  1964. 
Other  requests  shall  be  referred  to  U.S.  Naval 
Radiological  Defense  Lab . ,  Attn :  Weapons 
Capabilities  Branch,  Military  Evaluations 
Division,  San  Francisco,  CA  94135. 


_ AUTHORITY 

USNRDL  ltr  dtd  22  May  1967 


THIS  PAGE  IS  UNCLASSIFIED 


UNCLASSIFIED 

ad  4  4  1  9  8  3 

DEFENSE  DOCUMENTATION  CENTER 

FOR 

SCIENTIFIC  AND  TECHNICAL  INFORMATION 

CAMERON  STATION.  ALEXANDRIA.  VIRGINIA 


UNCLASSIFIED 


HOTICS:  When  governnent  or  other  drawings,  speci¬ 
fications  or  other  data  are  used  for  any  purpose 
other  than  In  connection  vlth  a  definitely  related 
governnent  proeurenent  operation,  the  U.  S. 
Qorernaent  thereby  Incurs  no  responsibility,  nor  any 
obligation  whatsoever;  and  the  fact  that  the  Oovern- 
aent  nay  have  fomilated,  furnished,  or  In  any  way 
supplied  the  said  drawings,  specifications,  or  other 
data  Is  not  to  be  regarded  by  implication  or  other¬ 
wise  as  In  any  Banner  licensing  the  holder  or  any 
other  person  or  corporation,  or  conveying  any  rights 
or  permission  to  manufacture,  use  or  sell  any 
patented  Invention  that  nay  in  any  way  be  related 
thereto. 


USNRDL-TR-741 
23  April  1964 


% 


C..D 

1  THE 


DEVELOPMENT  OF  A  WATER-SUBFACE-BURST  FALLOUT  MODEL: 
RISE  AND  EXPANSION  OF  THE  ATOMIC  CLOUD 


by 

I.  0.  Huebsch 


s 


not  to  be  released 

TJ.  S.  NAVAL  EADIOLOdftOQSAL 
DEFENSE  LABORATORY 


SAN  FRANCISCO  *  CALIFORNIA  *  94135 


*4 


4  19  83 


* 


* 


WEAPONS  CAPABILITIES  BRANCH 
S.  C.  Rainey,  Head 


MILITARY  EVALUATIONS  DIVISION 
S.  C.  Rainey,  Acting  Head 


ADMINISTRATIVE  INFORMATION 

This  work  was  sponsored  by  the  Defense  Atomic 
Support  Agency  under  NWER  Program  A-7a,  Subtask 
10.007.  It  is  further  described  as  Program  A -4, 
Problem  5,  in  the  USNRDL  Technical  Program  Summary 
for  Fiscal  Years  1964,  1965,  and  1966  (USNRDL- 
PR-46,  1  Oct  1963). 


DDC  AVAILABILITY  NOTICE 

Foreign  announcement  and  dissemination  of  this 
report  by  DDC  is  not  authorized. 


Eugene  R  Cooper 
Scientific  Director 


D.  C.  Campbell,  CART  USN 
Commanding  Officer  and  Director 


12ND  NRDL  PI  (9/63) 


1 


ABSTRACT 

A  water- surface-burst  fallout  model  is  being  developed  in  three 
parts:  (l)  cloud  rise  and  expansion,  (2)  particle  formation,  and  (3) 
particle  fall  through  the  atmosphere.  Part  (l)  is  discussed  in  detail. 
The  cloud  is  described  by  a  set  of  equations  for  rate  of  rise,  temper¬ 
ature  ,  water  content  and  other  parameters,  using  an  "entraining  parcel" 
method.  The  equations  allow  for  variations  in  atmospheric  temperature, 
pressure  and  humidity  as  well  as  in  initial  cloud  height,  temperature, 
and  air/water  energy  partition.  Temperature  and  volume  are  calculated 
by  extensions  of  standard  meteorological  equations.  The  key  assumptions 
in  the  cloud  model  are: 

1.  The  effective  rate  of  flow  of  environmental  air  into  the  cloud 
per  unit  surface  area  is  linearly  proportional  to  the  product  of  a 
characteristic  velocity,  and  the  ratio  of  cloud  density  to  ambient 
density. 

2.  Cloud  rise  is  retarded  by  both  entrainment  and  an  apparent 
eddy-viscous  force,  which  is  directly  proportional  to  the  same 
characteristic  velocity  and  inversely  proportional  to  the  density  ratio. 

3.  The  lost  kinetic  energy  of  rise  due  to  eddy  viscosity  remains 
in  the  cloud  as  turbulent  kinetic  energy,  diluted  by  entrained  air. 

4.  The  characteristic  velocity  involved  may  be  either  the  velocity 
corresponding  to  the  sum  of  the  kinetic  energies  of  cloud  rise  and  of 
turbulence,  or  the  greater  of  absolute  rate  of  cloud  rise  and  average 
velocity  of  turbulence.  Thus  entrainment  does  not  necessarily  end  when 
cloud  rise  ends,  but  continues  as  a  turbulent  diffusion. 

5*  The  time  required  for  the  cloud  to  accelerate  to  its  "terminal" 
velocity  from  rest  is  accounted  for  by  an  initial  virtual  mass  term. 

The  equations  have  been  programmed  for  machine  computation. 
Comparison  of  numerical  results  with  atomic-cloud  observations  indicates 
that  a  single  empirical  parameter,  the  entrainment  constant  \,  (which 
is  the  proportionality  constant  in  assumption  l)  can  also  be  used  as  a 
measure  of  the  rate  of  generation  of  turbulence.  Correct  predictions 
are  obtained  for  rate  of  rise,  cloud  size,  final  cloud  height  and  the 
late  horizontal  expansion  of  high  yield  clouds.  This  last  has  not  pre¬ 
viously  been  predicted  by  n  cloud  model.  Numerical  results  also 
indicate  that  for  low  yields,  the  effect  of  initial  water  content  of 
the  cloud  is  minor  compared  with  that  of  atmospheric  humidity  while  for 
high  yields  the  opposite  is  true. 

Finally,  the  particle-formation  and  particle-fall  portions  of 
the  model  development  are  summarized. 
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SUMMARY  PAGE 


The  Problem 

Land- surface  and  sea- water- surface  nuclear  explosions  produce 
different  kinds  of  fallout.  Water- sur face-hurst  fallout  particles 
consist  mainly  of  sea  salt  and  water.  There  is  no  fallout  model  for 
water-surface  bursts.  Until  now,  fallout  predictions  for  water- 
surf  ace-bursts  have  essentially  used  land- surface-burst  fallout  models. 

Findings 

NRDL  is  developing  a  water- surf ace -burst  fallout  model,  in  three 
parts.  Ibis  report  gives  a  detailed  discussion  of  Part  1:  Cloud  Rise, 
and  an  Outline  of  Work;  Part  2:  Particle  Formation,  and  Part  3s 
Particle  Fall  Through  the  Atmosphere.  The  cloud  rise  is  described  by 
a  set  of  equations  which  specify:  cloud  height,  size,  temperature, 
water  content,  etc . ,  as  functions  of  time  and  of  initial  and  * 

atmospheric  conditions.  The  equations  have  been  programed  for  machine 
computation  and  the  numerical  results  are  in  general  agreement  with 
weapon- test  observations  of  atomic  clouds. 
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SECTION  1 


INTRODUCTION 


1.1  FROBIEM  AND  OBJECTIVE 

Land- surface  and  sea-water-surface  nuclear  explosions*  produce 
different  kinds  of  fallout.  Sea-water- surface  burst  fallout 
particles  consist  mainly  of  sea  salt  and  water.  Their  water  contents, 
and  therefore,  their  masses,  change  during  their  fall.  Unlike  some 
land-surf  ace-burst  fallout  particles  (namely  unvaporized  soil 
particles)  they  do  not  exist  as  particles  before  the  burst. 

There  Is  no  fallout  model  for  water- surface  bursts.  Until  now, 
fallout  predictions  for  water- svrface  bursts  have  been  made  using 
land-surface  models,  with  at  most  minor  modifications.  (See  Sec.  1. 3) 

Therefore,  NRDL  Is  developing  a  water- surface-burst  fallout  model. 
This  report  gives  a  detailed  discussion  of  the  rise  and  expansion  of 
the  atomic  cloud  from  a  water-surface  burst,  and  an  outline  of  work 
on  the  other  parts  of  the  problem. 


- H - 

A  water- surface  burst  is  a  nuclear  explosion  centered  so  close  to 
the  water  surface  that  the  fireball  intersects  the  air-water  inter¬ 
face,  and  the  nuclear  cloud  initially  contains  both  air  and  vaporized 
water. 

The  terms  "atomic  cloud"  and  "nuclear  cloud"  are  used  interchangeably 
in  this  report. 
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1.2  BACKGROUND  AND  APPROACH 

Experimental  data  on  water- surf ace-buret  fallout  la  limited  to 
cloud-rise  histories  and  particle  collections  from  two  test  shots 
at  the  Paclfl  Proving  Grounds.  Both  shots  were  fired  on  barges  in 
relatively  shallow  water,  and  the  characteristics  of  the  fallout 
may  have  been  significantly  influenced  by  the  presence  of  barge 
and  bottom  material.  The  particle-collection  data  are  scattered, 
and  some  of  them  are  questionable.  These  data  are  insufficient 
for  constructing  an  empirical  fallout  model.  Therefore,  the  present 
model  is  being  developed  from  fundamental  physical  considerations 
using  the  available  experimental  data  as  guldeposts. 

The  development  has  been  divided  into  three  natural  inter¬ 
related  parte: 

1.  Tud  rise  and  expansion.  (The  main  subject  of  this 
report . ) 

2.  Fallout-particle  formation  in,  and  dispersion  from, 
the  clcmd. 

3>  Fallout-particle  fall  through  the  atmosphere,  subject 
to  meteorological  conditions. 

Preliminary  reports  on  the  three  parte  ere  being  published  separately. 
The  three  parts  will  then,  allowing  for  revision,  be  combined  in  the 
fallout  model  and  computer  program.  Refinement  and  changes  will  be 
considered  for  incorporation  into  the  model,  as  they  become  avail¬ 
able. 

Since  the  principal  objective  of  a  fallout  model  is  to  describe 
the  arrival  of  radioactivity  in  the  lower  atmosphere  and  at  the  sur¬ 
face,  Parts  1  and  2  of  the  development  exist  primarily  to  provide  a 
scientific  basis  for  this  description. 

This  is  a  research  fallout  model.  Numerical  computations  will 
show  which  variables  have  only  minor  effect  on  the  results.  For 
operational  use,  these  results  can  then  be  approximately  expressed 
in  simplified  equations,  or  in  tables,  graphs  or  nomograms. 


2 


v 


Another  summary  of  previous  cloud  models,  including  those  based 
on  the  work  of  Mechta1,  Sutton2*  3,  and  Taylor1*,  (see  also  the  dis¬ 
cussions  in  Refs.  5,  6,  7)  will  not  be  given.  Instead,  an  attempt 
will  be  made  to  satisfy  the  present  requirements  for  a  description 
of  cloud  rifle  (Sec  2.l). 

A  completely  analytic  approach  to  cloud  rise  seems  likely  to  bog 
down  in  hydrodynamics  and  to  result  in  delay  in  delivering  useful 
results.  Therefore,  this  report  uses  a  mixture  of  theoretical  analysis 
and  empiricism. 

1.3  POSSIBIE  INTERIM  WATER- SURFACE- BURST  FAILOUT  PREDICTION  METHODS 
The  following  interim  methods  have  been  used  or  considered  for 
prediction  of  water- surface-burst  (WSB)  fallout: 

1.  Use  a  land- surface-burst  (LSB)  model  such  ss  the  NRDL 

Q 

D-Model°  and  accept  the  inaccuracies  involved. 

2.  Make  small  adjustments  on  a  LSB  model  -  e.g.,  assume  a 
total  water- surface-burst  activity  equal  to  80$  of  that  from  an  LSB.  9 

3.  Use  the  D-Model  but  with  a  particle  size  distribution 
that  simulates  WSB  particle  arrival  times  and  locations  (at  sea  level), 
although  particle  mass  is  fixed.  THale  1.1  gives  typical  log-normal 
distributions  of  radioactivity  vb  particle  size,  used  with  the  D- Model 
for  Nevada  soil  and  coral  soil,  and  for  salt  particles. 

Table  1.1 

_ Log-Normal  Distributions  of  Radioactivity  vs  Particle  Size 

Mean  Particle 


Soil 

Diameter,  p 
(Microns) 

k>6±L_ 

Standard  Deviation 
of  Log  p 

Particle  Density 
(gm/cnH) 

Nevada 

113 

2.053 

0.732 

2.6 

Coral 

162 

2.209 

0.424 

2.6 

Salt  I 

122 

2.086 

0.177 

2.08 

Salt  II 

85*5 

1.932 

0.168 

2.08 
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In  computing  particle  falling  rates  using  the  "Salt"  distribution 
in  the  D- Model,  the  particles  are  taken  as  spheres,  while  Nevada  and 
Coral  particles  are  of  Irregular  shape  (cylinders  with  a  2*1  ratio  of 
length  to  width). 

1*.  Use  the  D-Model  rate  of  rise  and  expansion  but  consider 
variation  of  water  mass  (and  therefore,  of  falling  rate  at  a  given 
altitude)  of  slurry  particles,  using  the  computer  program10  developed 
for  Par low's  particle-fall  model.11  If  fallout  particles  consist  of 
water  and  salt,  we  must  not  only  specify  the  s&lt-mass  vs  activity 
distribution  but  also  initial  water  content.  (The  simplest  choice  is 
that  all  particles  are  initially  at  sea-water  concentration).  This 
question  is  considered  in  Sec.  4. 

5 .  Use  the  D-model,  but  replace  its  rate  of  rise  with  those 
calculated  for  a  water- surface  burst,  using  the  equations  of  the 
present  report. 

6.  Make  the  changes  in  the  D-model  described  in  both  4.  and 


5.  above. 


SECTION  2 


THE  RISE  OF  THE  ATOMIC  CLOUD: 

DISCUSSION  OF  PROBLEMS 

2.1  REQUIREMENTS  FOR  A  DESCRIPTION  OF  CLOUD  RISE 
The  Input  to  he  specified  consists  of: 

1.  Explosion  energy  and  energy  remaining  in  the  fireball. 

2.  Height  of  burst  (if  bursts  other  than  those  exactly  at 
sea  level  are  considered). 

3.  Atmospheric  temperature,  pressure,  and  humidity  as 
functions  of  altitude. 

ftie  desired  output  consists  of  the  following  cloud  properties  as 
functions  of  time: 

1.  Height. 

2.  Rate  of  rise. 

3.  Size  (vertical  and  horizontal  diameters). 

4.  Shape  (spherical,  ellipsoidal,  etc.). 

3>  Mass. 

6.  Volume. 

7.  Water  content  (vapor  and  liquid  or  solid). 

8.  Temperature  (either  mean  or  with  a  gradient  structure). 

9«  Turbulence  characteristics. 

The  cloud-history  output  is  needed  as  input  to  Part  2  of  the  fallout 
model:  particle  formation.  This  may  then  have  to  feed  back  to  Part  1, 
for  instance  if  the  formed  salt  particles  have  removed  a  large  part  of 
water  content  from  the  cloud  before  its  maximum  height  is  reached. 

The  effect  of  wind  on  cloud  rise  is  not  considered  in  this  report. 

It  may  be  considered  in  a  later  report  as  outlined  in  Sec.  4.1. 
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2.2  PARCEL  METHODS,  LOCAL  METHODS,  AND  SIMILARITY  RUIES 

2.2.1  Parcel  Methods  and  Local  Methods 

Parcel  methods  of  describing  the  rising  cloud  treat  the  cloud  as 
a  whole,  as  If  all  ports  of  the  cloud  had  the  same  properties  (temperature, 
density,  etc.).  Any  compensating  down  flow  of  air  outside  the  cloud  is 
neglected. 

Local  methods  give  equations  for  velocity  and  temperature  at  each 
point  in  the  flow  fields.  They  ore  mathematically  complicated  even 
when  the  field  has  axial  symmetry,  and  typically  describe  the  cloud 
by  a  set  of  partial  differential  equations.  (Such  equations,  of  course, 
refer  to  mean  quantities,  not  turbulent  fluctuations  at  each  point.) 

Not  only  is  it  simpler  to  use  on  "entraining  parcel"  to  represent  the 
cloud,  but  also  the  extreme  turbulence  prevailing  in  the  cloud  makes  it 
reasonable  that  the  effect  of  entrainment  is  immediately  averaged  over 
the  cloud.  Although  relative  velocities  of  parts  of  the  cloud  are  un¬ 
specified,  turbulence  parameters  can  still  be  used  to  represent  the 
net  effect  of  velocity  fluctuations. 

To  sura  up,  local  methods  require  numerical  solution  of  partial 
differential  equations,  unless  some  similarity  rules  (see  Sec.  2.2.2) 
can  be  devised.  Any  such  similarity  rule6  would  restrict  the  ability 
to  test  the  effect  of  variations  in  the  atmosphere  on  the  air/water 
composition  of  the  cloud,  on  specific  heat,  etc.  Such  testing  is  one 
of  the  objectives  of  this  "research  fallout  model."  Parcel  methods 
require  only  ordinary  differential  equations.  Therefore,  the  following 
discussion  (remainder  of  Sec.  2)  and  equations  (Sec.  3)  use  a  parcel 
method.  Use  of  a  parcel  method  does  not  keep  us  from  discussing 
entrainment  rules,  viscosity  analogies,  etc. 

2.2.2  Similarity  Rules 

If  two  sets  of  functions,  with  respect  to  both  position  and  time, 
of  velocity,  density,  tenperature,  etc.,  can  be  made  identical  by  an 
appropriate  choice  of  units  of  length,  time,  mass,  and  tenperature. 
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they  are  said  to  be  "similar"12  This  similarity  is  lacking  in  the 
rise  of  clouds  from  explosions  of  different  yields,  even  though  our 
parcel  method  replaces  the  continuous  flow  field  by  the  use  of  average 
cloud  and  ambient  values. 

1.  The  atmospheric  temperature -vs - altitude  gradient  changes 
at  the  tropopause  so  that  low-  and  high-yield  clouds  encounter 
different  conditions. 

2.  Even  if  the  atmosphere  has  a  constant  temperature 
gradient  (as  in  a  simple  model  of  the  troposphere)  the  large  temperature 
and  density  differences  between  cloud  and  ambient  air  violate  the 
usual  approximate  similarity  conditions.12  Also,  variation  of  specific 
heat  with  temperature  may  be  important. 

Therefore  no  attempt  is  made  in  this  report  to  scale  cloud  rise  from 
one  set  of  initial  and  atmospheric  conditions  to  another  (although  it  is 
assumed  that  the  11  characteristic  length”  involved  in  the  entrainment  and 
momentum- change  processes  discussed  in  Sec.  2.6  is  proportional  to  cloud 
radius).  Instead  the  cloud  rise  equations  (Sec.  3)  are  solved  numerically 
for  each  set  of  parameters.  Approximate  similarities  may  then  be  noted 
in  the  solutions. 

2.2*3  Potential  and  Virtual  Quantities 

The  use  of  potential  (adiabatic  equivalent) temperature  is  only 
possible  if  constant  specific  heat  is  assumed.  (Potential  temperature 
is  the  temperature  air  at  a  given  height  and  pressure  would  assume  if 
it  were  expanded  or  compressed  adiabatically  to  some  standard  pressure, 
say  sea- level  pressure.  If  the  effect  of  any  water  vapor  or  liquid  or 
solid  water  in  the  air  is  neglected,  the  corresponding  change  of 
temperature  with  height  is  called  a  dry-adiabatic  change.)  Also,  air 
entrainment  and  the  presence  of  water  vapor  complicate  the  use  of 
potential  temperature,  (if  the  presence  of  water  in  the  air  is  taken  into 
account,  the  adiabatic  change  of  temperature  with  height,  for  saturated 
air,  is  called  wet- adiabatic  or  saturated  adiabatic.  Because  of  the 
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release  of  latent  heat  of  condensation,  the  vet-adiabatic  rate  of  decrease 
of  temperature,  or  lapse  rate,  is  less  than  the  dry-adlabatic  lapse  rate. 
See  for  instance  Ref .  13 • ) 

It  is  convenient  to  express  some  cloud  properties  in  terms  of 
virtual  temperature  (the  temperature  of  dry  air  having  the  density 
of  the  actual  humid  air  at  the  given  temperature),  when  the  water- 
vapor  mass  fraction  of  the  air  under  consideration  is  6mall  (so  that, 
for  instance,  the  specific  heat  of  the  air  can  be  taken  as  that  of 
dry  air).  But  that  this  mass  fraction  is  small  is  precisely  what  ve 
cannot  assume  for  a  water- surface-burst  cloud.  We  still  need  the  actual 
temperature  of  the  cloud,  although  we  also  use  the  expression  for  the 
virtual  temper attire,  as  a  convenient  algebraic  abbreviation. 

2.3  INITIAL  CLOUD  CONDITIONS 

2.3*1  Thermodynamic  Conditions 

In  this  study,  consideration  of  the  fireball  begins  when  it  is  in 
approximate  pressure  equilibrium  with  the  atmosphere . ^5  (After  the 
"fireball”  starts  to  rise,  we  shall  call  it  the  "cloud.")  The 
temperature  at  this  time  (initial  temperature)  is  30O0°K  to  4000°K.  By 
not  going  back  to  higher  temperatures,  we  avoid  having  to  consider  air 
dissociation  and  loss  of  energy  by  thermal  radiation.  (The  cooling  of 
the  fireball  from  4000°K  to  3000°K  occurs  in  a  few  seconds  even  for  a 
megaton  shot. ) 

For  a  burst  centered  near  the  water  surface,  part  of  the  fireball 
gas  is  vaporized  sea  water  and  part  is  atmospheric  air.  Immediate  mixing 
of  the  air  and  salt- steam  volumes  is  assumed.  This  mixing  may  be  aided 
by  the  greater  buoyancy  (lower  molecular  weight)  of  water  vapor  rising 
into  the  hot  air.  Given  the  fireball  energy  and  temperature,  the  air- 
water  energy  partition,  and  the  specific  heats  of  air  and  water  as 
functions  of  temperature,  the  initial  cloud  radius  and  mass  can  be 
calculated  by  conservation  of  energy,  or  rather  enthalpy.  (Appendix  £.3, 
Initial  Conditions . )  It  is  assumed  that  the  dry  air  and  water  vapor 
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in  the  cloud  satisfy  the  perfect-gas  law,  but  that  their  specific  heats 
nay  vary  with  temperature. 

In  using  perfect-gas  theory,  we  are  not  claiming  the  fireball  was 
heated  reversibly;  of  course  it  was  not.  We  are  estimating  the  amount 
of  stored  energy.  To  use  an  analogy,  we  are  interested  not  in  the  fuel 
spilled  in  loading  a  rocket,  but  in  how  much  is  on  board. 

We  assume  the  initial  cloud  form  is  a  sphere.  For  a  water-surface 
burst  the  initial  conditions  could  center  the  "spherical"  fireball  at 
water  level,  or,  alternately,  at  a  height  equal  to  the  initial  cloud 
radius  so  that  the  fireball  Just  touches  the  surface ,  or,  more 
reasonably,  somewhere  between.  It  is  only  for  megaton  yields,  for 
which  the  initial  cloud  radius  is  a  kilometer  or  more,  that  this  var¬ 
iation  in  initial  height  of  cloud  center  is  large  enough  to  have  a 
significant  effect  on  final  cloud  height  and  size. 

2.3*2  Air-Water  Energy  Partition 

Hie  calculated  initial  cloud  size  is  not  very  sensitive  to  the 
assumed  air- water  energy  partition.  An  all-air  cloud  has  about  10$ 
greater  radius  than  one  with  a  50-50  energy  partition,  because  the 
greater  specific  heat  of  water  and  the  large  latent  heat  of  evaporation 
overcompensate  for  the  smaller  molecular  weight  of  water  as  compared  to 
that  of  air,  and  result  in  a  smaller  fireball. 

One  might  expect  that  a  check  on  the  air-water  partition  could  be 
made  by  observing  the  initial  size  and  rate  of  rise  of  the  cloud  from  a 
water- surface  burst  as  compared  with  a  land- surface  burst  or  a  low  air 
burst  of  the  same  yield.  But  calculations  with  the  equations  in  Sec.  3 
indicate  that  for  a  megaton  surface  burst,  the  differences  in  cloud 
height  and  radius  of  a  50-50  burst  as  compared  with  those  of  a  100$- air 
burst  near  the  surface,  are  too  small  to  be  observed  on  photographs,  if 
the  differences  are  due  solely  to  the  energy  partition.  Therefore,  any 
observed  differences  must  be  due  to  other  factors.  (Different  initial 
teiqperatures  are  one  possibility. ) 


9 


2.4  EFFECT  OF  WATER  CONTENT  ON  THE  CLOUD 

2.4.1  Effect  of  Water  Content  on  Cloud  Rise 

The  Importance  of  the  Bea- water  content  of  the  fireball  is  due  tos 

1.  Hie  accompanying  Bea  Balt  which  condenses  and  coagulates 
to  form  fallout  particles  (second  part  of  this  problem). 

2.  The  possible  additional  thrust  the  cloud  receives  relatively 
late  in  its  rise  from  release  of  latent  heat  of  condensation,  and  of 
freezing.  Once  the  rising  cloud  has  cooled  to  below  condensation 
temperature,  the  temperature  -  height  curve  (neglecting  entrainment) 

is  wet-,  not  dry-,  adiabatic  (Sec.  2.2.3)>  so  that  the  cloud  then  cools, 
and  loses  buoyancy,  more  slowly.  Energy  release  by  condensation  appears 
to  be  significant  only  for  low-yield  clouds  (most  of  whose  water  comes 
from  the  atmosphere).  Calculations  indicate  that  megaton  clouds  reach 
the  stratosphere  before  cooling  to  the  boiling  point,  so  the  extra  thrust 
may  come  suxpdsingly  late.  But  the  water  fraction  of  the  cloud  is  then 
so  small  as  to  make  the  thrust  insignificant,  unless  the  cloud  has  a 
hot  core  or  ring,  which  is  not  considered  here. 

Hie  influence  of  the  salt  content  of  sea  water  ( 3$  by  mass )  on 
cloud  density,  temperature,  and  water-vapor  pressure  is  neglected. 

Atmospheric  humidity  also  influences  cloud  rise.  Computations 
with  the  cloud-rise  equations  given  later  in  this  report  indicate  for 
a  1  KT  cloud,  entrained  atmospheric  water  vapor  is  far  more  important 
than  initial  water  content,  while  for  a  1000  KT  (l  IfT)  cloud,  the 
opposite  is  true.  Hie  reason  for  this  difference  is  that  the  cloud 
surface-to-volume  ratio  decreases  as  yield  and  cloud  size  increase,  so 
that  entrainment,  which  is  proportional  to  this  ratio  (Sec.  2.6.1), 
becomes  leBB  influential. 

Computations  indicate  that  water-vapor  condensation  begins  below 
the  freezing  point  (273°K)  for  megaton  clouds,  so  that  the  latent  heat 
released  may  include  that  from  the  llquid-to-ice  transition.  A  choice 
as  to  which  latent  heat  to  use  is  complicated  by  the  presence  of  salt 
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and  by  turbulence,  both  of  which  influence  the  freezing  point.  One 
must  ask  the  question:  does  the  15#  additional  latent  heat  of  freezing 
influence  cloud  rise  much?  This  will  be  tested  by  computation. 

Possibly  the  radar  image  of  a  water  surface  burst  cloud  would 
look  "denser"  than  that  of  a  non-water  buist.  This  is  an  interesting 
subject  for  future  investigation. 

2.4.2  Effect  of  Water  Loss  on  Cloud  Buoyancy 

In  the  later  stages  of  cloud  rise,  much  of  the  water  content  is  in 
liquid  or  solid  form.  The  question  arises:  does  all  liquid  and  solid 
water  remain  in  the  cloud,  or  is  some  of  this  "spent  fuel"  Jettisoned, 
lightening  the  cloud's  load? 

If  the  water- and- salt  particles  coagulate  to  a  size  with  a  sensible 
falling  rate  while  the  cloud  is  still  rising  (this  will  be  discussed  in 
detail  in  a  later  report),  then  there  is  a  gain  in  cloud  buoyancy  by  loss 
of  these  particles.  Hie  maximum  possible  influence  of  this  effect  is 
measured  by  the  ratio  of  liquid  and  solid  water  mass  to  dry  air  mass 
in  the  cloud.  A  small  value  of  the  water  percentage  of  the  fireball 
energy  partition,  together  with  a  large  entrainment  rate  and  a  low- 
humidity  atmosphere,  would  minimize  the  possible  gain  in  buoyancy.  If 
this  fallout  from  the  cloud  takes  place  after  the  end  of  cloud  rise,  it 
should  have  little  effect  on  cloud  buoyancy  because: 

1.  The  cloud  has  gained  so  much  air  by  entrainment  that  the 
water  contribution  to  total  density  is  small. 

2.  The  cloud  no  longer  moves  as  a  unit.  See  Ref.  l6,  (Sec.  13, 
Settling  of  Clouds)  for  a  criterion  for  this  mass- subsidence  problem. 

2.4.3  Effect  of  Water  Loss  on  the  Cloud's  "Heat  Engine" 

Can  an  MT  cloud,  rising  past  its  equilibrium  level  in  the  strato¬ 
sphere  against  a  positive  ambient  temperature  gradient,  have  lost  so 
much  of  the  condensed  (and  coagulated)  water  by  fallout,  or  have  entrained 
so  much  dry  air,  that  it  is  then  cooling  dry-adiabatlcally?  Calculations 
indicate  this  can  occur,  but  since  dry  nd  wet  adiabatics  hardly  differ, 
at  these  low  temperatures,  the  effect  -n  cloud  motion  i6  negligible. 
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2.5  INFLUENCE  OF  VARIATION  OF  AMBIENT  CONDITIONS 

The  volume  and  mass  of  a  rising  cloud  change  by  two  processes: 

1.  Adiabatic  expansion. 

2.  Entrainment  of  arnbl  'nt  air. 

How  sensitive  to  variations  In  meteorological  conditions  are 
these  processes?  The  first  proc  ss  depends  on  atmospheric  pressure, 
which  does  not  vary  much  (at  a  g^v  i  altitude).  For  high-yield  clouds, 
most  of  the  second  process  would  occur  at  greater  altitudes  where  time 
variations  of  meteorological  conditions  are  smaller  than  In  the  lower 
troposphere.  But  as  mentioned  In  Sec.  2.4.1  low-yield  clouds,  on  the 
other  hand,  are  sensitive  to  ambient  humidity. 

While  theories  of  atomic  cloud  behavior  have  largely  been  suggested 
by  studies  of  cumulus  cloud  behavior,  there  ~e  limitations  on  the  use 
of  such  meteorological  theories  to  describe  atomic  clouds.  Cumulus 
clouds  have  been  assumed  to  consist  of  a  series  of  "thermals"  (buoyant 
convective  elements).  The  atomic  cloud  Is  a  single  thermal.  Some  orders 
of  magnitude  are  quite  different  for  atomic  clouds  and  cumulus  clouds. 

The  Initial  temperature  difference  between  a  cumulus  cloud  .and  its 
environment  Is  about  1°C  or  less.  This  Is  enough  to  move  the  cloud  up, 
and  even  bounce  It  off  the  tropopause  and  back  down  several  times, 
according  to  some  theoretical  calculations . ^  Nuclear  clouds  have 
tremendous  Initial  temperature  differences,  much  greater  speeds  of 
rise  and,  for  megaton  bursts,  much  greater  sizes  than  cumulus  clouds. 

Thus,  atmospheric  temperature  conditions  are  much  less  Important 
In  the  early  stages  of  cloud  rise  for  nuclear  clouds  than  for  cumulus 
clouds.  If  the  cloud  temperature  is  1000°K,  It  does  not  Immediately 
matter  whether  ambient  temperature  Is  290°K  or  300°K.  The  assumed 
Initial  rate  of  rise  of  the  cloud  makes  a  big  difference  In  the  sub¬ 
sequent  rise  of  a  "thermal"  on  the  meteorological  scaled  but  has 
negligible  Influence  on  an  nuclear  cloud.  Hence,  we  can  assume  that  a 
nuclear  cloud  (or  fireball)  Is  still  at  rest  when  It  has  cooled  to  our 
Initial  temperature  (about  3000°K).  Also,  unlike  the  case  of  cumulus 
cloud  rise,  ambient  turbulence  can  be  neglected,  and  surely  so  can 
"erosion"  (evaporation  at  the  surface  of  the  cloud)  except  for  very  low 
yields  or  after  stabilization. 
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2.6  A  MODEL  OF  ENTRAINMENT  AND  MOMENTUM  CHANGE 

Several  approaches  to  cloud  rise  were  considered  in  the  course  of 
this  study.  These  approaches  included: 

Static  methods: 

Thermal  equilibrium 
Potential-kinetic  energy  balance 

Dynamic  (momentum- equation)  methods: 

Potential  flow 
Drag  coefficient 

External  or  internal  momentum  exchange. 

The  various  approaches  are  discussed  in  Appendix  B.  The  internal- 
momentum-  loss  rule  and  the  closely  related  entrainment  rule  that  were 
adopted  are  discussed  in  this  section. 

2.6.1  Entrainment 

It  was  suggested  by  Taylor  that  the  rate  of  increase  of  cloud  mass 
due  to  entrainment  of  air  is  proportional  to  surface  area  times  absolute 
rate  of  cloud  rise: 


dm 

dt 


S  \ 


(2.1) 


where  m  is  cloud  mass 
t  is  time 
S  is  surface  area 

A  is  a  dimensionless  constant  of  proportionality 
dz 

u  *  is  rate  of  rise,  where 


z  is  cloud  (center)  height 
J?e  is  ambient  (environment)  density 
This  equation  is  obtained  by  assuming  that,  averaging  over  the  cloud 
surface,  8,  the  ambient  air,  of  density j?e  ,  flows  "onto"  (or  into) 
the  cloud  at  a  rate  proportional  (factor  \)  to  the  absolute  value  of 
rate  of  rise,  |u|  .  (The  absolute  value |u|  used  here  and  later  avoids 
negative  entrainment  or  drag,  if  the  cloud  sinks  down  after  overshooting 
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its  equilibrium  height.)  For  a  spherical  cloud  of  radius  r,  this 
leads  to  dr  =  X  |dz|  if  the  change  in  cloud  density  is  neglected. 
Empirical  values  of  X  for  nuclear  and  other  clouds  range  from  0.07 
(Ref.  7)  to  0.30  (Ref.  5)  and  0.3^  (Ref.  18).  Theoretical  and 
experimental  work  on  entrainment  has  considered  mostly  cumulus 
clouds  and  turbulent  Jets,  at  about  the  same  density  and  temperature 
as  ambient  air.  In  actuality,  the  rising  cloud  is  initially  about 
ten  times  as  hot  and  one  tenth  as  dense,  a6  the  ambient  air.  We 
propose,  therefore,  that  the  entrainment  rate  is  also  proportional 
to  the  "relative  inertia"  of  cloud  and  ambient  air,  i.e.,  the 
density  ratio j Ifet  (where y-is  the  cloud  density)  so  that. 


dm 

dt 


mX  I  v 


(2.2) 


whore 

V  **  cloud  volume 

v  is  a  characteristic  velocity  given  by: 


V  =  V  u2  +  EE. 

where  Ejr  is  the  turbulent  energy  density  (Sec.  2.6.2).  This 
entrainment  rate  agrees  qualitatively  with  the  remark  in  Ref.  18 
that  "in  the  early  stages  of  its  motion  the  atomic  bomb  cloud 
behaves  more  as  a  bubble  of  air  in  water  in  which  case  there  is  no 
entrainment  at  all."  Also,  weighting  the  entrainment  rate  by  the 
density  ratio  avoids  the  difficulty  encountered  by  the  authors  of 
Ref.  7.  So  as  not  to  cool  the  cloud  too  fast  they  had  to  make  X 
"surprisingly  small."  Then  the  choice  of  initial  radius  required 
to  make  the  final  cloud  size  correct  was  found  to  be  "surprisingly 
large"  and  was  described  as  an  "artificial  radius." 

Any  Justification  for  introducing  the  density  ratio  is  necessarily 
heuristic,  but  an  argument  for  this  approach  can  be  made  in  mechanical 
terms.  Suppose  that  a  lump  of  cloud  gas,  with  velocity  u,  collides  with 


a  lump  of  stationary  ambient  air  and  Imparts  a  velocity  to  it,  and 
suppose  that  the  corresponding  momentum  and  kinetic  energy  are  con¬ 
served  in  a  collision.  Momentum  is  not  conserved  in  the  overall 
cloud- atmosphere  interaction,  since  turbulence  is  generated  at  the 
expense  of  directed  motion,  but  ve  can  imagine  the  Interaction  as 
consisting  of  separate  momentum-transfer  and  momentum-dissipation 
processes.  Suppose  now  that  a  lump  of  ambient  air  must  gain  a 
certain  threshold  velocity  ku  to  be  entrained,  (where  k  is  a  con¬ 
stant),  and  that  the  lump  of  cloud  gas  remains  in  the  cloud.  We 
calculate  the  volume  of  ambient  air  Ve  to  which  a  unit  volume  of 
cloud  can  transfer  this  threshold  velocity.  (Actually  the  cloud 
lump  may  make  many  collisions,  but  we  idealize  the  process  to  a 
single  collision. )  The  equations  for  conservation  of  momentum  and 


kinetic  energy  are,  respectively, 

u  = pn'  +fe  Ve  ku  (2.3) 

\  ve  k2  u2  (2.4) 

where  u'  is  the  velocity  of  the  cloud  lump  after  collision. 

The  solution  is  Vg  =(^>^e)(2/k  -  l)  (2.5) 

u'  =  u  (k  -  l)  (2.6) 

(This  solution  has  no  physical  meaning  for  k>2.) 


Thus,  the  volume  (and  consequently  the  mass)  of  ambient  air, 
that  a  lump  of  cloud  of  unit  volume  can  entrain  with  a  given 
threshold  velocity,  is  proportional  to  the  density  ratio  y  /pe  . 

We  could  also  say  that  the  distance  that  the  lump  travels  before 
losing  its  momentum  is  proportional  top /pe,  and  call  this  distance 
the  "characteristic  length."  This  length  could  then  be  considered 
the  distance  between  the  inside  and  the  outside  of  the  cloud. 
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In  the  particular  case,  k  ■  1,  the  threshold  velocity  is  u 
itself,  the  entraining  cloud  lump  is  left  with  zero  velocity,  and 
the  entrained  volume  Jj>Q  times  the  entraining  volume,  so  that 
the  entrained  mass  equals  the  entraining  mass.  Now  suppose  that  in 
time  dt  a  volume  of  cloud  gas  S  x|v|dt  "flows  out"  of  the  cloud. 

(This  is  the  Taylor  hypothesis^  applied  to  the  entraining  air  instead 
of  the  entrained  air.)  Then  the  entrained  mass  is 


dm  = S  X  v  dt. 


so  that 


dm  _  S 
dt  V 


mX  Iv 


dt 


as  was  suggested  above. 


(2-7) 


(2.2) 


More  formally,  we  can  estimate  the  entrainment  rate  by 
regarding  momentum  as  a  diffusible  quantity.  In  a  turbulent  flow, 
however,  ordinary  molecular  diffusivity  is  of  minor  importance  com¬ 
pared  with  the  diffusive  effect  of  turbulence.  As  has  been  done  by 
many  investigators,  we  suppose  the  diffusive  effect  of  turbulence 
can  be  expressed  as  an  "eddy  diffusivity"  having  the  same  dimensions 
as  molecular  diffusivity,  although  the  eddy  diffusivity  is  a  pro¬ 
perty  of  the  flow  conditions,  not  of  the  fluid.  Just  as  molecular 
diffusivity  is  the  product  of  the  speed  and  mean  free  path  of  a 
molecule,  times  a  constant,  we  may  take  the  eddy  diffusivity  to  be 
proportional  by  a  factor  X  ,  to  the  product  of  the  speed,  v,  of  a  lump 
of  cloud,  and  a  length,  /£,.  We  call  v  and  the  "characteristic  speed" 
and  "characteristic  length,"  respectively. 

Suppose  the  characteristic  speed,  v  to  be  the  speed  corresponding 
to  the  total  kinetic  energy  of  a  lump  of  cloud.  This  energy  consists 
of  the  kinetic  energy  of  rise  and  the  kinetic  energy  of  turbulence,  E^, 

so  that  _ 

v  *  V  u2  +  Ek  (2.8) 

where  the  positive  value  of  the  square  root  is  taken. 
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Suppose  the  characteristic  length,  JZ>,  to  he  proportional  to  cloud 
(vertical)  radius,  and  that  this  length  is  the  distance  from  inside 
to  outside  the  cloud.  Then,  assuming  that  the  rate  of  change  of 
momentum  per  unit  volume  with  distance  along  the  length  has  the 
average  valued  v/  Qm  the  total  outward  flow  of  momentum  in  time  dt 
is  S (kvjfr)  (^v/£)  dt,  (where  X  v£/  is  the  eddy  diffusivity),  so 
that  (dividing  the  momentum  flow  in  dt  by  the  momentum  per  unit  volume 
^7v)  the  corresponding  change  in  volume  is  SXvdt,  as  assumed  above  in 
Eq>  (2.2).  The  boundary  of  the  cloud  moves  outward,  as  a  result  of  this 
increase  in  volume.  As  outlined  in  Sec.  2.6.3,  entrainment  continues 
after  the  end  of  cloud  rise,  when  there  is  no  longer  any  net  total 
cloud  momentum,  and  now  v  consists  entirely  of  the  average  velocity 
of  turbulence, V  2E^ 

The  following  analogy  to  the  entrainment  process  was  proposed  by 
C.  F.  Ksanda.  Consider  the  rising  cloud*  of  radius  r,  as  a  piston  of 
the  same  radius.  Let  the  cloud  rise  during  an  Interval  dt.  Hie  re¬ 
sulting  environment  process  consists  of  two  steps. 

(l)  Neglecting  expansion,  an  element  of  cross  sectional 
area  dA  at  the  top  of  the  cloud  rises  a  distance  dz, 
and  sweeps  out  a  volume  dAdz^. 


(2)  Suppose  that  the  ambient  air  in  this  volume  is  engulfed 
by  the  cloud.  Its  mass  is  j?edAdz-^.  If  the  entrained  mass 
is  small  compared  go  the  cloud  mass  m,  it  now  expands  from 
dAdz^  to  cLAdz^ ),  using  the  gas  law  at  constant  pressure. 

Hie  total  Increase  in  cloud  height  in  time  dt  is  thus  not 
dz^  but  dz  =  (^g  ^)dz^. 


17 


Representing  the  top  surface  of  the  cloud  by  its  horizontal  projection 
A y  the  total  entrained  mass  is  dm  = ^Adri-j>  and  we  can  write 


dm  A 
dt  “  m  V 


u 


(2.10) 


For  a  sphere,  A  -  s/4  where  S  is  the  surface  area.  This  analogy 
thus  predicts  \  =  .25,  which  is  actually  in  good  agreement  with 
observation.  The  analogy  does  not  of  course  apply  in  the  late 
stages  of  entrainment  when  the  driving  foice  is  not  cloud  rise  but 
internal  turbulence,  but  at  that  time  cloud  and  ambient  densities 
are  nearly  equal. 

2.6*2  Momentum  Loss  and  Turbulent  Kinetic  Energy  Gain 

To  determine  momentum  loss,  we  draw  an  analogy  between  turbulence 
and  molecular  motion,  as  we  did  for  entrainment.  In  this  case,  how¬ 
ever,  we  speak  of  eddy  viscosity  instead  of  eddy  diffusivity. 

Suppose  that  the  rising  cloud  is  retarded  by  the  stationary 
ambient  air,  as  if  the  latter ,  like  the  cloud,  were  an  "eddy- viscous" 
fluid,  with  the  same  eddy  viscosity  \  v  JL .  Then  the  force  on  a  unit 
area  of  cloud  surface  is  proportional  to  y^\vj£(u/£/),  where  the  actual 
velocity  gradient  is  replaced  by  an  average  value  u/L.  The  total  force 
is  then  proportional  to  Sj?v^(u/^  )  and  the  force  per  unit  cloud  mass 
to  (S/v)(pe^ )  \  v  u. 

It  may  be  argued  that  not  the  total  surface  area  but  the  projected 
area  normal  to  the  direction  of  rise  should  be  used.  For  a  sphere,  the 
ratios  of  these  quantities  to  V  are  respectively  3/r  and  3/4r.  Further¬ 
more,  we  do  not  know  whether  eddy  diffusivity  and  eddy  viscos¬ 
ity  are  the  same,  or  merely  proportional.  Therefore,  for 
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generality,  we  replace  S/V  by  a  constant  times  l/r,  combine  the  con¬ 
stant  with  \  and  write  the  product  as  2kg,  so  that 

force  per  unit  cloud  mass  „  2k0  \  v  u  (2.11) 

/  “ 

Computations  indicate  that  when  this  form  is  used  the  empirically 
correct  value  of  2k2  is  very  nearly  the  same  jl/  appearing  in  the 
entrainment  equation. 

If  the  retarding  force  on  the  cloud  were  an  ordinary  viscous 
force,  the  kinetic  energy  loss  due  to  this  force  would  appear  as 
Increased  heat,  i.e.  random  motion  of  molecules.  Since  it  is  an 
eddy- viscous  force,  the  lost  energy  appears  as  an  increase  in  turbulent 
kinetic  energy,  i.e.  random  motion  of  lumps  of  fluid,  which  are  sup¬ 
posedly  mixed  into  the  cloud.  Thus,  eddy  viscosity  does  not  reduce 
the  total  cloud  kinetic  energy,  but  converts  some  of  the  energy  of 
rise  to  turbulent  energy. 

Then,  neglecting  entrainment,  the  rate  of  Increase  of  turbulent 
kinetic  energy  per  unit  mass  is  equal  to  the  rate  of  decrease  of  kinetic 
energy  of  rise  due  to  eddy  viscosity, 


dEfc 

dt 


k2 


\  v  u2 
/  r 


(2.12) 


We  regarded  the  eddy  viscous  momentum  loss  as  occurring  on  the 
outside  of  the  boundary  of  the  cloud,  in  fluid  of  density  j^.  The  same 
final  result  could  be  obtained  by  regarding  the  momentum  loss  as 
occurring  on  the  inside  of  the  boundary,  in  fluid  of  density P,  and  with 
a  different  velocity  gradient  given  by  u/ (  Jl This  would  agree  with 
the  comment  in  Sec.  2.6.1  that  the  distance  that  a  lump  of  cloud  travels 
while  transferring  momentum  is  proportional  to  p/pe  •  This  suggests  that 
if  we  follow  the  original  Taylor  entrainment  rule,  Eq.  (2.1)  in  Sec.  2.6.1, 
we  should  then  omit  the  facto TP &  in  the  eddy- viscous  force. 
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That  is,  the  present  model  gives  more  drag  and  less  entrainment  by  the 
factors  fa  If*  and  respectively,  than  the  "original"  rule.  These 

arguments  are,  however,  suggestive  rather  than  convincing,  anu  equally 
plausible  ones  might  be  made  to  obtain  other  results. 

Qualitatively,  the  net  effect  on  cloud  rise  of  introducing  the 
density-ratio  hypothesis  iS! 

1.  The  cloud  cools  more  slowly  since  it  entrains  relatively 
less  low-altitude  air  and  more  high-altitude  air  whose  potential 
temperature  (Sec.  2.2.3)  is  higher. 


2.  The  early  rate  of  mass  increase  is  less. 

3.  There  is  relatively  less  drag  by  entrainment  and  relatively 
more  by  eddy  viscosity,  so  that  more  turbulence  is  generated. 


We  can  now  equate  the  rate  of  change  of  cloud  momentum?  per  unit 
mass,  to  the  sum  of  the  forces  acting  or.  tfc-2  cloud. 


du  /e  -P 
dt  ~  f 


g  -  2k2 


& 


|v|u 


/ 


1  dm 
m  dt 


u 


(2.13) 


Here  the  three  terms  on  the  right  are  respectively  due  to  buoyancy 
(i.e.  Archimedian  force),  drag  and  entrain nent .  In  the  momentum 
equation  in  Sec.  3*1  the  r  in  the  denominator  of  the  drag  term  is 
replaced  by  a  generalized  characteristic  leigth  X,  which  assumes 
different  values  for  spherical  and  spheroidal  clouds. 

If  the  total  drag  force  on  a  spherical  cloud  is  represented  by 
a  drag  coefficient  Cp,  then 


CD 


F 


D 


it  r2(^u2/2) 


(2.1k) 


» 
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so  that  the  drag  force  per  unit  cloud  mass  is  given  by 


Vm  -  °D  (2.15) 

which  has  exactly  the  form  of  the  drag  term  in  the  momentum 
equation,  Eq.  (2.13). 

Finally,  we  can  allow  for  the  time  for  the  cloud  to  accelerate 
from  rest  to  terminal  velocity  by  treating  the  initial  motion  of  the 
(spherical)  cloud  as  a  potential  flow,  which  sets  in  motion  a  volume 
of  ambient  air  equal  to  one-half  the  initial  displaced  volume.  All 
entrained  fluid  becomes  turbulent,  so  that  the  initial  cloud  is  pre¬ 
sently  "dissolved  in  turbulence."^  Therefore  we  multiply  the  entire 
right  side  of  the  momentum  equation  by  ,  where  m'  is  a  virtual 

mass  equal  to  the  mass  of  one  half  the  initial  displaced  volume.  This 
virtual  mass  factor  is  derived  in  App.  B.2.1. 

2.6.3  Additional  Effects  of  Turbulence 

When  cloud  rise  ends,  no  more  kinetic  energy  can  be  transferred  to 
turbulence.  The  turbulence  drives  the  apparent  horizontal  expansion 
of  the  cloud  which  continues  (for  high-yield  clouds)  after  the  end  of 
the  rise,  vertically  stabilized  by  the  stratospheric  temperature 
inversion  (See  Sec.  2.8).  In  our  parcel  method,  using  the  character¬ 
istic  velocity: 

v  =  V  u2  +  2Ejj 

in  the  entrainment  equation,  Eq.  (2.2),  this  late  expansion 
corresponds  to  the  case  u  *  0,  v  =  V  2Ek 

If  instead  we  chose  as  characteristic  velocity, 

v  =  max  (|u|,\/2Ek  ),  then  the  idealized  picture 

would  be  that  as  cloud  rise  slows  down,  the  diffusive  spread  of  tur¬ 
bulence  (turbulence-induced  mixing)  overtakes  the  radial  entrainment 
flow  induced  by  the  mean  flow. 
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During  cloud  rise  and  expansion,  the  turbulence  influences 
coagulation  of  fallout  particles  in  the  cloud. -**9  Thus,  retardation  of 
cloud  rise,  late  horizontal  expansion  of  cloud,  and  particle  growth 
give  us  three  check  points  for  the  constant  kg>  (Sec.  2.6.2),  which 
governs  the  rate  of  transfer  of  kinetic  energy  from  directed  motion 
to  turbulence.  As  noted  in  Sec.  2.6.2,  computation  indicates  that 
2kg  is  related  to  the  entrainment  parameter  X,  so  that  X  is  the 
only  independent  empirical  quantity  in  the  description  of  the  cloud. 

It  is  assumed  that  only  a  small  fraction  of  the  turbulent  energy  is 
degraded  to  heat  by  viscosity  during  cloud  rise,  i.e.  relaxation  time 
is  large  compared  to  cloud  rise  time.  Computations  of  the  late  hor¬ 
izontal  expansion  of  megaton  clouds  confirm  this  hypothesis. 

2.6.4  Compressibility  Effects  on  Cloud  Rise 

For  yields  over  15  MT,  the  present  equations  (Sec.  3),  as  well 
as  potential  flow  methods,  predict  cloud-rise  speeds  that  may  be  too 

high,  if  the  drag  coefficient  or  eddy- viscosity  coefficient  used  is  the 
same  as  that  used  for  lower  yields.  For  a  100  MT  cloud,  speeds  as 
high  as  sonic  speed,  at  ambient  temperature,  are  predicted. 

We  have  assumed  that  the  "characteristic  length"  increases 
linearly  with  cloud  radius.  But  at  some  large  radius  and  speed  of  rise, 
it  would  seem  that  a  small  part  of  the  boundary  no  longer  knows  what 
the  cloud  diameter  is,  i.e.  has  time  to  adjust  itself  to  an  increase  in 
diameter.  That  is,  the  time  for  information  to  pass  half  way  around 
the  periphery  is  not  small  compared  with  the  time  for  the  cloud  to  rise 

a  significant  distance,  say  one  radius,  i.e.  —  is  at  least  as  large 

r  u  1  c 

as  of  say  —  ,  where  c  is  the  speed  of  sound. 

u’  c  jt 

Compressible  flow  theory  suggests  multiplying  the  eddy-viscous 
drag  by  a  factor  —  — — ,  where  M  is  the  Mach  number,  M  «  u/c. 

\J  1  -  M2 

(A  more  sophisticated  correction  would  permit  M  =  1.)  This  suggestion 
is  heuristic,  and  based  on  calculations  for  slender  bodies  at  small 
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angles  of  attack.  Furthermore,  even  M  ««  2/3  only  produces  a  35$ 

Increase  In  drag,  using  this  factor. 

The  local  speed  of  sound  near  the  boundary  of  the  cloud  Is  higher 
than  In  the  ambient  air,  because  of  the  high  cloud  temperature  so  that 
the  local  Mach  number  is  smaller,  and  the  "correction"  less  effective. 
But  also,  the  local  flow  speed  near  the  cloud  equator  may  be  greater 
than  the  net  upward  speed,  increasing  the  Mach  number.  These  comments 
are  offered  only  as  notes  for  future  investigation. 

2.7  VORTICITY  AND  CIRCULATION 

Hie  details  of  vortex  structure  and  cloud  circulation  are  not  con¬ 
sidered  in  this  study.  We  are  not  interested  in  vortex  structure  but 
in  the  total  turbulent  energy  arising  from  eddy-viscous  momentum  loss 
during  cloud  rise.  We  put  this  energy  in  a  "black  box"  and  let  it 
contribute  to  the  characteristic  velocity  given  by  Eq.  (2.8)  and  thus 
to  the  entrainment  process.  Just  what  vortices  are  in  the  black  box 
does  not  matter.  In  fact,  the  box  contains  the  turbulent  spectrum. 

We  do  not  use  a  vortex  model  such  as  that  of  Levine. ^  In  his 
model  the  bubble  exchanges  mass  with  the  environment.  The  narrowness 
of  the  wake  of  atomic  clouds  Indicates  very  little  such  exchange  takes 
place.  And  even  this  narrow  wake  may  just  be  environmental  air  lifted 
above  the  condensation  level  without  mixing.  For  the  nuclear  cloud 
case  there  is  no  inactive  portion  of  the  cloud  through  which  the  bubble 
is  rising,  since  the  whole  cloud  is  the  bubble.  Unlike  a  cumulus  cloud, 
it  is  not  true  that  "the  turbulence  is  very  likely  confined  to  a  thin 
layer  the  thickness  of  which  is  small  compared  with  that  of  the  bubble." 
Furthermore,  Levine's  model  uses  an  exchange  factor  (See  App.  B.2.3). 

In  the  nuclear  cloud,  all  the  retarding  ambient  air  is  entrained,  (except 
possibly  the  "initial  virtual  mass"  in  our  model)  so  there  is  no  loss 
of  total  kinetic  energy  by  exchange,  but  only  transfer  from  mean  motion 
to  turbulence. 
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2.8  SPHERICAL  AND  NON- SPHERICAL  CLOUDS 

We  assume  that  the  initial  cloud  form  is  a  sphere.  If  and  when  the 
top  of  an  atomic  cloud  reaches  or  passes  the  tropopause  it  is  observed 
to  flatten  (by  horizontal  expansion,  not  vertical  contraction)  into  an 
oblate  spheroid.  We  represent  this  flattening  by  postulating  that 
vertical  expansion  ~ase s  when  the  cloud  center  reaches  the  tropopause, 
so  that  the  cloud  becomes  an  oblate  spheroid.  The  entrainment  rule  is 
unchanged.  (Actually,  even  for  clouds  not  passing  the  tropopause 
horizontal,  diameter  grows  somewhat  more  than  vertical  diameter.)  What 
are  the  mechanics  of  this  shape  change? 


1.  Effects  of  the  finite  size  of  the  cloud.  The  top  and  base 
of  the  cloud  encounter  different  conditions.  The  larger  the  cloud,  the 
less  valid  it  is  to  approximate  ambient  conditions  by  those  at  altitude 
of  the  cloud  center.  Near  the  tropopause,  the  atmosphere  has  a  non¬ 
linear  lapse  rate.  When  the  cloud  reaches  the  tropopause,  its  top 
loses  buoyancy  more  rapidly  than  its  bottom  which  is  still  below  the 
tropopause.  Picturing  the  cloud  as  an  elastic  body  meeting  increased 
resistance  to  its  rise,  we  can  say  that  to  conserve  momentum  it  can 
only  move  outward.  This  suggests  that  the  rate  of  distortion  is  pro¬ 
portional  to  the  difference  in  buoyancy  between  its  top  and  bottom;  or 
better,  to  the  rate  of  change  of  this  difference  with  altitude.  That 


is,  a  sharp  increase  in  relative  stability  in  the  atmosphere,  or 


roughly, 


>  0, as  at  the  tropopause,  flattens  the  cloud.  (Here  Te 


is  ambient  temperature . )  Speculatively, 


<0  would  produce  a 


vertical  elongation  of  the  cloud  into  a  prolate  spheroid.  This 

PI 

speculation  is  supported  by  a  recent  meteorological  cloud  rise  study. 

A  canputation  was  made  for  a  conditionally  unstable  atmosphere,  and 
"unlike  the  spherical  shape  envisaged  in  bubble  theories,  the  cloud... 
develops  into  a  tall  and  slender  current."  So  we  should  say  that 
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distortion  rate  is  proportional  to  rate  of  change  of  stability.  But 
for  a  nuclear  cloud  only  a  very  sharp  stability  change,  such  as  at 
the  tropopause,  can  affect  cloud  form. 

2.  Turbulent  diffusion.  For  a  high-yield  cloud,  the 
accumulated  turbulent  kinetic  energy  induces  turbulent  (horizontal) 
diffusion  of  cloud  particles  with  the  vertical  thickness  fixed  by  the 
stably  stratified  atmosphere.  The  rate  of  horizontal  "expansion" 
depends  on  the  turbulent  energy  density  and  thus  on  the  rate  of  eddy- 
viscous  momentum  loss  during  rise.  (See  Sec.  2.6.2.)  Low  yield  clouds 
rise  more  slowly,  generating  less  turbulence,  and  their  energy  density 
decreases  rapidly  by  dilution  by  entrained  air  so  that  such  clouds 
hardly  expand  after  the  end  of  rise. 
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SECTION  3 


A  SET  OF  EQUATIONS  FOR  CLOUD  RISE 


The  following  set  of  equations  for  cloud  rise  uses  an  "entraining 
parcel"  method  (Sec.  2.2)  and  the  entrainment  and  eddy- viscosity  con¬ 
cepts  discussed  in  Sec.  2.6.  The  numbered  equations  are  those  forming 
the  set,  which  have  been  solved  numerically  (Sec.  3«1*0  •  The  unnumbered 
equations  are  used  only  in  deriving  the  numbered  equations .  The  alter¬ 
nate  (a)  equations  correspond  to  the  hypothesis  that  the  "characteristic 
length,"  and  therefore  the  rates  of  momentum  change  and  entrainment,  do 
not  depend  on  the  density  ratio  (Sec.  2.6).  The  set  includes  "dry"  (d) 
and  "wet"  (W)  equations  for  the  water  vapor  content,  temperature,  volume, 
and  liquid  or  solid  water  content  of  the  unsaturated  and  saturated  cloud, 

respectively.  The  momentum,  water  vapor  and  temperature  equations  are 

17  22 

suggested  by  those  used  in  parcel-method  cumulus  cloud  calculations.  ’ 

3.1  MOMENTUM 

As  in  Eq.  (2.13),  we  equate  the  rate  of  change  of  momentum  to  the 
sum  of  buoyancy  and  eddy- viscous  forces, 

d  2k2  v  u 

dt(mu)  =  v  (Pe  “  p)g - —m 

where  t  is  time 

m  is  cloud  mass 
u  is  rate  of  cloud  rise 
V  is  cloud  volume  =  m/p 
p  is  cloud  density 
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Pe  is  ambient  air  density 
g  is  acceleration  due  to  gravity 
kg  is  constant 

,  vhere  is  the  turbulent  energy  density 
^  is  a  "characteristic  length"  (Sec.  2.6) 

Using  the  perfect  gas  law  for  air  and  water  vapor  since  the  pressure 
inside  and  outside  the  cloud  is  the  same,  then  with  appropriate  weights. 


where  T  and  Tg  are  respectively  cloud  and  ambient  temperature 
q(x)  may  be  shown  to  equal 

x  and  xg  are  respectively,  cloud  and  ambient  "mixing  ratios" 

(ratios  of  water- vapor  mass  to  dry  air  mass  per  unit  volume) 
w  is  the  ratio  of  liquid  and  solid  water  mass  to  dry  air 
mass  in  the  cloud 

c  «  18/29,  the  ratio  of  the  molecular  weight  of  water  vapor  to 
that  of  dry  air 


1  +  x 
1  +  x  +  w 


- =  P  is  the  ratio  of  cloud  gas  density  to  total  cloud 
density.  (This  term  allows  for  the  liquid  and  solid  water 


content  of  the  cloud,  if  any.  The  effect  of  this  content 
is  usually  small. ) 


We  use  the  symbols:  T  =  Tq(x)  and  T  =  T  q(x  ). 

#  #  e  e  e  ^ 

T  and  T  are  thus  the  virtual  temperatures  (Sec.  2.2.3)  and  T  P  is 
e 

the  (cloud)  virtual  temperature  allowing  for  the  contribution  of  liquid 
and  solid  water  to  the  total  cloud  density. 


Taking  the  characteristic  length  proportional,  to  the  density 
ratio  so  that  =  £  Pe/Pe>  where  (  depends  only  on  cloud  size,  then 
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the  momentum-balance  equation  may  be  rewritten  as 


* 

r 

T 

— =  0  -  1 
Te* 

g  - 

2k2  vT* 


„  .  1  dm 
*  ra  dt  u 


(3-1) 


Alternatively,  if  is  independent  of  p/pe  so  that  =  ^, 


"  -x 
du  T 

dt  m 

1e 


P  -  is 


[*2kg  v 

l8-nr  + 
_ 


(3.IA) 


Since  in  potential- flow  theory,  energy  and  momentum-balance  statements 

require  that  the  fluid  mass  accelerated  from  rest  include  a  virtual 

mass  equal  to  one  half  the  initial  displaced  mass,  then  the  entire 

right  side  of  each  of  Eqs.  (3«l)  and  (3*1A)  must  be  multiplied  by 

- 2-jg - —  ,  where  the  subscript  o  indicates  the  initial  value  of 

m  +  mo  T072Te0* 

each  quantity.  (See  App.  B.2.1.) 

3.2  HEIGHT 


The  height,  z,  of  the  center  of  the  cloud  is  given  by 

dz 

dt  ~  u 

3.3  WATER  VAPOR 
3.3.1  Dry 


(3.2) 


During  the  "dry"  (unsaturated)  period  no  water  is  lost  by  condensa¬ 
tion*  In  time  dt,  a  mass  dm  of  ambient  air  is  entrained,  so  that,  at 
time  t  +  dt,  the  new  "mixing  ratio"  will  be  given  by 


x(t  +  dt)  = 


x  ,  ,  e 
m  *= — ; —  +  dm  *= — ; — 

1  +  X  1  +  > 


1  +  X  1  +  X 


28 


The  specific  heat  of  entrained  air  is  taken  as  that  of  dry  air, 

c  (T).  For  a  small,  interval  dt,  the  sum  of  the  following  two  terms  is 
pa 

zero  (by  conservation  of  energy): 
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Heat  used  to  warm  entrained  air  dm: 


and 


cpa(T)  dT' 


dm 


Heat  absorbed  in  adiabatic  expansion  dp: 


m  [cp(T)  dT  -  RT*  &] 

where  cp(T)  is  the  weighted  mean  of  the  specific  heats  of  air  and 
water  vapor, 


(T)  +  cpw(T) 
1  +  x 


The  expansion  term  becomes,  using  the  hydrostatic  equation  and  the  gas 
law, 

* 

m  [cp(T)  dT  +  gdzj 

e 

Equating  the  sum  of  the  two  changes  to  zero, 


dT 

dt 


A  v<t)  * 


1  dm 
m  dt 


(3.4D) 


The  two  terms  on  the  right  give  the  effects  on  temperature  of  adiabatic 
expansion  and  entrainment,  respectively. 


3.4.2  Wet 


Since  the  temperature 


of  the  saturated  cloud  is  at  most  373°K, 
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specific  heats  are  taken  as  independent  of  temperature .  As  before, 

is  the  weighted  mean  of  c  _  and  c  .  For  a  small  interval  dt,  the  heat 

pa  pw  2 

released  by  condensation:  -mldx/(l  +  x)  ,  is  equal  to  the  sum  of 
Heat  used  to  saturate  and  warm  entrained  air: 


L(x  -  x  )  dm  +  c  „  (T-T  )  dm 

e  Pa  e 


and 


Heat  absorbed  in  expansion  dp: 


m 


c  dl«  R  T* 

P  a  p  : 


This  heat  balance  neglects  the  small  heat  release  by  cooling  of  liquid 
or  solid  water,  -mwcdX,  where  c  -  specific  heat  of  liquid  or  solid  water* 
Performing  the  indicated  balance,  as  for  the  "dry"  heat  equation  and 

g 

also  using  Eq.  (3.3W),  and  dropping  small  terms  of  order  x  ,  and  taking 
cpa  ^  cp» 


dT 

dt 


1 


(x-xj  ~  +  (T  -  Tjiig 


1  + 


HZ 


e  ;  m 

i 


(3.4W) 


The  two  terms  on  the  right  (corresponding  to  the  two  terns  in  the 
numerator )  are  respectively  the  ("wet")  adiabatic  and  entrainment  terms. 

3*5  VOLUME 

In  a  nmaii  height  dz,  cloud  volume  changes  an  amount  dV,  due  to: 


1.  Adiabatic  expansion  at  constant  mass,  dV)m,  and 
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2.  Maes  entrainment  at  constant  pressure,  dV)p,  so  that 
dV  =  dv)m  +  dV)p 

To  calculate  each  dV,  the  gas  law,  p  =  2  r  q(x)  T,  is  used  in  the 
differential  form: 

dp  dm  t  dT  ,  dq(x)  dV 

p  ■  -  *  T  *  -JtJT  •  T 

This  assumes  the  non-gas  mass  fraction  is  small. 


3.5.1  Die 


c  (T)  is  the  temperature-dependent  weighted  mean  of  c  and  c  . 
p  pa  pv 

For  an  adiabatic  volume  change  at  constant  mass,  dm  and  dq  are  zero 

y 

--  ”,  (where  7  is  the  ratio  of  specific  heats  of  the  cloud 
gas),  so  that 


and  — ^ 

p 


dV)m  dT)m  dT)ffl 
V  *  T(y-l)  "  T 

Using  the  adiabatic  term  in  Eq.  (3*4d), 


R  qfx)  “ 


dV), 


1  1  ' 
m  *  *  1  R.  q(x)  "  cpj 


1 sated 


For  constant-pressure  entrainment, 
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dT)^  is  given  by  the  entrainment  terra  in  Eq.  ( 3 • ) ,  and 


dq 


Adding  the  two  volume  changes , 


l/€-l 

(1  ♦  x)2 


dx 


«-» 


e  q(x)  1 

1  i\  | 

W  1 

=Pj 

,  T 

/  c  (T)  dT 

J  Pa' 

T 

,  e  ll  dm  /„ 

1 - -= - im  at  <3-5D) 


u 


Tc 


p 

-1  dx 


i(x)  (l  +  xf  at 

3.5.2  Wet 

Using  the  adiabatic  term  in  Eq.  (3*4w), 

f  1  *  FT  1 

a  i 


=  ^  rbr-r 

T  \  aMA  p 


1  + 


L2€x 


cRT 
P  a 


2  : 


Using  the  entrainment  term  in  EH.  (3*^W), 

r  -  i 


dv)p  =  v  1 1  - 


(x-xe)  —  +  (T-Te) 


1  + 


L2€X 

cRT2) 
P  a 


dm 

m 


where  dq  is  neglected. 
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Adding  the  tvo  volume  changes, 


* 

X 

1 

)  h  .  (T  -  T  ) 
c  +  e 

P 

1  dm 

T2*  \ 

m  dt 

T 

c  R  T 

P  a  /  J 

3.6  LIQUID  AND  SOLID  WATER  CONTENT 


(3.5W) 


Liquid  and  solid  vater  content  appear  only  in  the  wet  (w)  equations. 
If  the  cloud  expands  or  warms  up  so  much  that  its  vapor  pressure  drops 
below  saturation,  then  the  dry  equations  again  apply.  The  "wet" 
equations  do  not  know  about  this.  But  if  the  ratio  of  liquid  and  solid 
water  mass  to  dry  air  mass,  w,  drops  to  zero  in  Eq.  (3.6W),  a  return 
to  the  "dry"  equations  is  required. 

3.6.1  Dqt 

Let  w  be  the  ratio  of  liquid  and  solid  water  mass  to  dry  air  mass, 
w  =  Then, 


w  =  0 


(3-6d) 


3.6.2  Vet 

The  liquid  and  solid  water  mass  can  increase  by 
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1.  Excess  of  the  mixing  ratio  of  entrained  air  over  that  of 
the  saturated  cloud,  and 

2.  Condensation  of  vapor  already  in  the  cloud,  so  that 

% i  V*  an  .  to 

dt  1  +  x  dt  "  a  dt 


E|y  definition  of  w. 


m. 


a.  . 

iat  ■ '  v,jr  +  _5t 


V  +  X  -  x_ 


1  +  X 


dm 


dx 


dt  ma  dt 


Substituting  m 


m 


1  +  x  +  v 


dv 

dt 


1  ;  1  +  x 

“  S  1 1  +  x 
>  ef 


V  +  X  -  X. 


1  dm  dx 
m  dt  ”  dt 


3.7  TURBUIEMT  KINETIC  ENERGY  DENSITY 


(3.6m) 


The  energy  transfer  rate,  or  the  kinetic  energy  loss  corresponding 
to  the  kg  (momentum  loss)  term  in  Eq.  (3*1),  is  the  rate  of  generation 
of  turbulent  energy  per  unit  mass  (Sec.  2.6.2,  Eq.  (2.12)).  The  gain 
of  turbulent  energy  per  unit  mass  is  diluted  by  the  increase  of  mass, 
so  that  the  turbulent  energy  density  E^  is  given  by 


dt 


=  k. 


* 

T 

¥ 

Te 


2 

u  v  p  1  dm 
K  m  dt 


(3.7) 
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Alternatively,  it  f  is  independent  of  p/p  , 

e 


dt 


k2u2v 

T~ 


F  1  dm 
k  m  dt 


(3-7A) 


3,8  MASS 

As  discussed  in  Sec.  2.6.1, 


1  dm 
m  dt 


(3.8) 


where  S  is  cloud  surface  area.  This  is  Eq.  (2.2)  in  Sec.  2.6.1. 
Alternatively,  if  the  entrainment  rate  is  independent  of  density  ratio, 


1  dm 
m  dt 


S  _T _ 

V  * 
T 

e 


S  \  v 


(3.8A) 


This  follows  from  Eq.  (2.1)  in  Sec.  2.6.1. 
3.9  CLOUD  FORM:  VOLUME 


For  a  spherical  cloud 


(3-9) 


where  r  is  the  radius.  If  a  certain  altitude  z^,  (T  for  tropopause) 
is  designated  such  that  the  vertical  radius  remains  equal  to  r(z^)  = 
after  2  =  then  for  the  resulting  oblate  spheroid 


r  =  I3SL. 
h  ; 


\  1/2 


(3.9E) 


36 


vhsre  r.  ie  the  horizontal  radius  (major  semi-axis)  of  the  spheroid. 


3.10  CLOUD  FORM:  EFFECTIVE  SURFACE  AREA 


For  a  spherical  cloud 


S  =  4nrc 


(3.10) 


FOr  an  oblate  spheroid  cloud 


r.  1 

2  .  rT  1  +  e 


S  *  2*  r.  +  n  -  -  In 
n  v 


(3.10B) 


where  e  is  the  eccentricity, 


2  2 
e-  ~  rT 

rh 


3.11  CHARACTERISTIC  IENCTH 

For  the  present  we  disregard  compressibility  effects  (2.6.4)  and 
take  as  characteristic  length,  % 


(3.11) 


for  the  sphere,  and 


i  -  rT 


(3.HE) 


for  the  spheroid. 


3.12  THE  EFFECT  OF  FREEZING 


To  allow  for  the  effect  of  freezing,  let  L  be  the  latent  heat  of 
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sublimation,  when  T  is  less  than  some  temperature  Tp.  This  neglects 
the  liquid- to- solid  heat  release  if  condensation  occurs  both  above 
and  below  Tp  and  the  difference  in  vapor  pressure  over  ice  and  water 
below  Tp. 

3.13  COMMENTS  ON  THE  SET  OF  EQUATIONS 

The  set  of  equations  to  be  solved  consists  of  the  eight  differential 
equations  (3»l)  to  (3*8),  with  the  three  algebraic  equations  (3*9)  to 
(3.11)  (or  (3.9E)  to  (3.UE)  for  a  spheroid)  substituted  in  them. 

The  initial  conditions,  defining  equations  for  parameters  and  values 
of  constants  are  given  in  Appendix  C. 

The  total  number  of  equations  can  be  reduced  by  substitution.  But 
the  present  extended  form  is  more  easily  revised  end  more  suitable  for 
machine  computation. 


Certain  equations  can  be  divided  by  u  =  ^  to  give  a  "quasi-static" 
subset  with  z  as  the  independent  variable.  Under  such  restrictive 
assumptions  as  constant  specific  heat,  zero  water  content  and  the 
characteristic  length  rule,  f ^  independent  of  density  ratio,  as  in 
Eq.  (3.1A),  an  analytic  solution  for  r(z)  can  be  obtained.  These 
algebraic  curiosities  axe  not  discussed  further  in  this  report. 

3.14  NUMERICAL  SOLUTION  OF  THE  EQUATIONS 

The  eight  differential  equations  (3*1)  to  (3*8)  with  the  defining 
equations  (3.9)  to  (3«H)  have  been  programmed  for  machine  computation, 
using  a  Runge-Kutta  method.  The  program  permits  the  alternative 
density  ratio  assumption  (use  of  the  "A"  equations)  and  certain  other 
options.  Appendix  E  gives  the  FORTRAN  program,  using  the  characteristic 
velocity  v  =  max  (|u|,  V  2  Ej^)  as  discussed  in  Sec.  2.6. 3« 
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Tables  3*1  and  3*2  are  sample  computer  outputs  for  two  sets  of 

2q 

initial  conditions  and  parameter  values.  A  modified  tropical  atmosphere 
was  used.  Appendix  D  gives  a  glossary  of  printout  symbols.  All  quanti¬ 
ties  are  in  mks  units,  except  explosion  energy  W  in  kilotons. 

Table  3.1  is  a  sample  output  for  a  5  MT  explosion  with  l/3  of  the 
energy  remaining  in  the  initial  cloud,  centered  to  Just  intersect  sea 
level.  Note  the  late  horizontal  expansion  of  the  cloud  (column  R). 

Computations  were  also  made  for  the  same  initial  conditions  and 
parameters  as  those  for  Table  3*1>  tut  using  the  alternate  equations 
(3>1A),  (3.7A),  and  (3*8a),  omitting  the  initial- virtual-mass  factor, 
and  with  k^.  =  0,  i.e.  no  eddy-viscous  force  and  thus  no  turbulence  effects. 

^  4 

This  is  essentially  the  original  Taylor  hypothesis  except  for  the 
sphere- spheroid  transition  at  the  tropopause.  It  was  found  that  this 
hypothesis  did  not  account  for  the  late  horizontal  expansion  of  the 
cloud,  and  predicted  much  larger  vertical  oscillations  of  the  cloud, 
after  the  cloud  reached  maximum  height . 

Table  3*2  is  the  output  for  a  low  yield  (20  KT)  shot.  It  shows 
the  much  lower  turbulent  energy  density  (column  EK)  of  low  yield  clouds. 
Note  also  that  the  cloud  does  not  reach  a  maximum  height  at  about  360 
sec.,  but  continues  rising  slowly.  This  is  an  effect  of  the  humid 
tropical  atmosphere  ("conditional  instability")  and  agrees  with  Pacific 
test  observations. 

The  numerical  results  are  in  general  agreement  with  observations 
of  atomic  cloud  radius,  rate  of  rise  and  (for  megaton  clouds)  late 
horizontal  expansion,  when  kg  is  about  0.10  to  0.15  and  X  about  0.20 
to  0.30  (and  when  the  initial  virtual  mass  factor  is  included  in  the 
momentum  equation).  This  indicates  that  (as  suggested  in  Sec.  2.6.2) 

2kg  is  closely  related  to  X,  so  that  2kg  may  be  replaced  by  X. 
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Then  the  cloud  model  vould  have  only  one  empirical  parameter,  X,  vhoee 

4  18 

magnitude  is  found  to  be  that  given  in  earlier  ctudies.  ' 

Preliminary  testing  of  the  effect  of  atmospheric  conditions  on 
cloud  rise  indicates  that  for  low  yields,  the  initial  water  content 
of  the  cloud  is  unimportant  compared  with  atmospheric  humidity.  For 
very  high  (megaton)  yields  the  opposite  is  true,  but  the  water  content 
has  only  minor  effect. 
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SECTION  4 


OUTLINE  OF  FUTURE  WORK  ON  THE 
WATER- SURFACE- BURST  FALLOUT  MODEL 

4.1  THE  NUCLEAR  CLOUD 

As  suggested  in  Sec.  2.6.4,  a  compressibility  correction  for  the 
rise  of  very  high  yield  clouds  may  be  introduced. 

The  effect  of  wind  shear  on  cloud  rise  may  be  considered  as  follows: 
As  the  ratio  of  the  wind  shear  between  top  and  bottom  of  cloud  to  the 
rate  of  rise  increases,  the  cloud  is  stretched,  increasing  the  surface- 
to-volume  ratio  and  therefore  the  entrainment  rate,  and  damping  the 
cloud  motion. 

The  cloud-rise  characteristics  predicted  by  the  equations  using 
different  parameter  values,  will  be  compared  with  nuclear  cloud 
observations. 

4.2  PARTICLE  FORMATION 

The  atomic  cloud  from  a  sea- water- surface  burst  contains  a  large 
amount  of  water,  a  much  smaller  amoirrt  of  sea  salt,  and  a  very  small 
amount  of  metallic  material . 

As  the  cloud  cools,  the  metal, salt  and  water  condense  successively. 
The  condensed  metal  particles  are  the  nuclei  for  particle  growth.  The 
metal- salt-water  particles  are  called  slurry  particles  or  slurry  drop¬ 
lets.^ 

We  will  assume  that  radioactivity  is  distributed  uniformly  per  unit 
sea  salt  mass.  Given  an  estimate  of  metal  and  salt  masses  per  unit 
cloud  volume,  the  initial  particle  masses  can  be  taken  arbitrarily  small 
as  the  rate  of  coagulation  by  Brownian  movement  is  sensitive  to  the  mass 
density,  not  number  density  for  small  particles.  (That  is,  it  does  not 
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matter  if  we  suppose  the  mass  to  consist  of  one  particle  of  diameter  10“^ 
cm  or  a  thousand  particles  of  diameter  10  cm  per  unit  volume,  because 
of  the  rapid  growth  from  10"?  to  10“^  cm). 

As  the  particles  reach  about  1  micron  (l0“^cm)  diameter  the 
coagulative  effect  of  Brownian  movement  decreases  and  coagulation  by 
turbulence^  must  influence  further  particle  growth.  Part  of  this 
growth  occurs  during  the  dry  stage  of  cloud  rise  and  part  during  the 
wet  (including  frozen)  stage. 

The  temperature,  liquid  and  solid  mass  per  unit  cloud  volume,  and 
intensity  of  turbulence  all  vary  with  time.  Hence,  approximations  will 
be  needed  in  analyzing  the  coagulation  process,  in  order  to  use  existing 
steady- state  theories. 

It  may  be  that  large  slurry  particles  are  formed  only  from  megaton 
water- surf ace- bursts.  This  follows  if  the  rate  of  generation  of 
turbulence  is  proportional  to  the  cube  of  rate  of  rise,  as  in  Eq.  (3*6), 
Sec.  3*6.  This  suggests  that  low-yield  water-surface  bursts  would 
produce  relatively  little  fallout,  since  the  requisite  coagulative 
influence  would  be  lacking. 

4.3  FALLOUT  OF  SLURRY  PARTICLES 

Until  now^®>  ^  slurry-droplet  mass-altitude-time  histories  have 
been  calculated  with  the  particles  starting  their  fall  either  dry  or 
wet  (at  sea  water  concentration)  from  a  reasonable  stabilized  cloud 
height  (say  60,000  ft  for  a  megaton  burst)  and  from  various  lesser 
altitudes . 

When  calculations  for  "initially  wet"  (i.e.  at  sea-water  salt 
concentration)  particles  were  compared  with  weapon  test  data^'  ^  it 
was  found  that  particles  with  large  salt  mass  came  down  too  soon  and 
those  with  small  salt  mass  came  down  too  late,  with  some  crossover  point 
at  which  calculated  and  observed  times  of  arrival  agreed. 

A  tentative  solution  to  the  discrepancy  is  as  follows:  in  the 
earlier  part  of  cloud  rise,  particles  are  being  formed  dry.  As  they 
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grow,  some  of  them  escape  from  the  cloud  by  diffusion  or  gravity.  The 
largest  particles  have  the  greatest  probability  of  escape.  Since  they 
escape  dry,  they  are  Initially  lighter,  and  reach  the  ground  later  than 
wet  particles  with  the  same  salt  mass.  The  remaining  "salt  solution" 

In  the  cloud  Is  diluted  by  the  loss  of  salt  to  below  sea-water  salt 
concentration.  (Entrained  atmospheric  moisture  also  contributes  to 
this  effect.)  Thus,  the  smaller  particles  which  stayed  In  the  cloud 
larger  come  out  wetter.  Then  they  are  larger  and  heavier,  per  unit  salt 
mass,  and  reach  the  ground  earlier  than  the  "Initially  wet"  particles 
of  Refs.  10  and  11,  and  In  better  agreement  with  weapon  test  data. 

Par low10  has  offered  a  similar  argument. 


SECTION  5 


RESULTS  AND  CONCLUSIONS 


5.1  RESULTS 

A  model  of  the  rise  end  expansion  of  the  atomic  cloud  from  a  water- 
surface  hurst  has  been  developed.  A  set  of  differential  equations 
describes  the  changes  with  time  in  cloud  size,  shape,  altitude, 
temperature,  water  content  and  turbulent  energy  density  as  functions 
of  explosion  energy,  initial  air-water  energy  partition  and 
atmospheric  conditions.  The  model  is  applicable  to  air  bursts  as 
well  as  water  surface  bursts.  It  uses  turbulence  concepts  which 
acoount  for  the  late,  horizontal  expansion  of  megaton  clouds,  and  can 
also  be  applied  to  the  formation  of  fallout  particles.  The  equations 
have  been  programmed  for  machine  computation.  Numerical  results  in¬ 
dicate  that  with  appropriate  choice  of  parameters,  the  cloud  model 
agrees  with  observations  of  nuclear  explosion  clouds.  This  agreement 
will  be  the  subject  of  a  separate  (Classified)  report. 

5.2  CONCLUSIONS 

The  cloud  model  gives  an  adequate  phenomenological  description  of 
the  rise  and  expansion  of  a  nuclear  cloud.  This  description  can  pro¬ 
vide  thermodynamic  input  data  for  development  of  the  second  part  of 
the  water- surf ace  buret  fallout  model:  particle  formation.  It  can  also 
specify  the  initial  altitude  of  particles  for  Part  3  of  the  model: 
particle  fall  through  the  atmosphere.  Part  1  of  the  three-part 
development  of  the  fallout  model  is  therefore  accomplished.  Since 

Part  3  of  the  model  has  already  been  developed,  only  Part  2  remains  to 
be  developed,  followed  by  combination  of  the  three  parts  into  the  final 
model. 
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APPENDIX  A 


SYMBOLS  USED  IN  THE  REPORT 

A.l  A  NOTE  ON  NOTATION 

This  report  uBeB  hydrodynamics,  thermodynamics  and  meteorology. 
These  fields  use  the  Bame  symbols  for  different  quantities; 
consequently,  any  notation  uBed  muBt  violate  some  usage.  For  example, 
in  meteorology  x  and  w  are  UBed  for  ratiOB  of  vapor-  and  liquid-water 
mass  to  dry  air  maBB,  respectively.  But  in  hydrodynamics  the  velocity 
components  u,  v,  w  correspond  to  the  coordinates  x,  y,  z.  Since  z 
is  the  usual  symbol  for  the  vertical  coordinate,  as  in  dp  =  -peg  dz, 
inconsistency  cannot  be  avoided. 

A. 2  SYMBOLS  USED  IN  THE  PRESENT  REPORT 

a  horizontal  Bemi  axis  of  cloud 
b  vertical  Bemi  axis  of  cloud 

Op  drag  coefficient 

c  specific  heat  of  liquid  or  solid  water;  speed  of  sound 

Cp  specific  heat  of  gas  at  constant  pressure 

Cy  specific  heat  of  gas  at  constant  volume 

Eft  turbulent  kinetic  energy  per  unit  mass 

eB  saturation  vapor  pressure  of  water 

f  fraction  of  explosion  energy,  W,  contained  in  fireball  at 
start  of  rise 

g  acceleration  of  gravity 

relative  humidity 


6o 


empirical  constant  (in  eddy  viscosity) 
latent  heat  of  evaporation 


characteristic  length 
molecular  weight;  Mach  number 
mass  of  cloud 

pressure 
1  +  x/c 
1  +  x 

universal  gas  constant  * 

p 

gas  constant  of  air  =  — — 

Ma 

1  +  x/c 

weighted  mean  value  of  gas  constant  =  Ra  “ — 

radius  of  cloud 
surface  area  of  cloud 
temperature 

Tq(x ),  i.e.  virtual  temperature 


vertical  velocity 
volume  of  cloud 

characteristic  velocity,  v  =  yu^+2  or  max(|u|,  -/?P'K) 

total  explosion  energy  (kilotons) 
liquid  and  solid  water  mass  per  unit  dry  air  mass 
mixing  ratio  (water  vapor  mass  per  unit  dry  air  mass) 
vertical  coordinate 

lapse  rate  of  atmosphere  =  rate  of  decrease  of  temperature 
with  height 

ratio  of  gas  density  to  total  density  of  cloud  =  ,1  +.  x _ 

1  +  x  +  w 


ratio  of  specific  heats  Cp/cy 

ratio  of  molecular  weights  of  water  and  air  =  18/29 

empirical  constant  (in  entrainment  rate) 

density 


6l 


fraction  of  fireball  energy  used  to  heat  air 

SUBSCRIPTS 

a  air  (dry  air) 

e  ambient  (environment)  conditions 

h  horizontal  (radius  of  cloud) 

o  initial  condition 

v  water  or  water  vapor 

wv  water  vapor 

M £■  liquid  and  solid  water  (i.e.  water  and  ice) 
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APPENDIX  B 


CLOUD  HEIGHT,  CLOUD  VELOCITY,  AND  ENTRAINMENT: 

POSSIBLE  APPROACHES 

A  set  of  cloud- rise  equations  must  contain  a  momentum  equation  if 
cloud  properties  are  to  be  calculated  as  functions  of  time.  Without  a 
momentum  equation,  some  of  these  properties,  such  as  size  and  tempera¬ 
ture  can  be  determined  as  functions  of  height  only,  and  dynamical 
properties  such  os  turbulence  characteristics  cannot  be  given  at  all. 

By  a  static  or  quasi-static  description  of  cloud  rise  we  mean  one 
lacking  a  momentum  equation. 

B.l  STATIC  AND  QUASI- STATIC  METHODS 

Static  and  quasi-static  methods  give  a  maximum  cloud  helght--l.e. 
tell  at  what  altitude,  but  not  when,  to  stop  cloud  rise. 

B.1.1  Static  Equilibrium 

In  the  static-equilibrium  approach,  cloud  rise  stops  when  cloud 
temperature  equals  ambient  temperature.  Actually  the  cloud  may  have 
enough  momentum  to  rise  far  beyond  thl6  height,  meanwhile  entraining 
more  air. 

B.l. 2  Quasi-Static  Approach  or  Energy  Balance 

While  the  cloud  Is  hotter  or  lighter  than  ambient,  it  loses  poten¬ 
tial  energy  and  gains  kinetic  energy  as  it  rises;  when  it  is  cooler  or 
heavier  the  opposite  energy  change  takes  place  (for  a  saturated  cloud, 
latent-heat  energy  release  is  also  considered). 

Using  this  approach,  Ref.  2k  calculates  temperature  in  a  non¬ 
entraining  "giant  thunderstorm"  cloud  (observed  top,  70,000  ft)  by 
balancing  the  two  areas  between  the  ambient  and  the  cloud  (wet  adiabatic) 
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temperature  vs  height  curves,  namely  the  areas  above  and  below  the 
intersection  of  the  curves#  The  cloud  picks  up  kinetic  energy  while 
it  is  warmer,  and  loses  it  (above  the  intersection)  while  it  is  cooler. 
This  gives  maximum  height.  For  an  (entraining)  nuclear  cloud,  the 
areas  would  have  to  be  weighted  for  variation  of  cloud  mass  with  height. 

Weighting  the  buoyancy  vs  height  function,  <^(z),  according  to 
cloud  mass  m,  the  cloud  stops  at  the  height  z^  for  which 

m(z)  £p(z)dz=  0  (B.l) 

Here  »  (pe  -  p)/pe,  where  pe  and  p  are,  respectively  ambient  and 
cloud  densities.  Now  the  cloud  is  denser  than  the  environment.  To 
find  an  equilibrium  height,  one  brings  the  cloud  down  again  to  z2, 
for  which 

z2 

m  Zsp  dz  =  0  (B.2) 

1 

and  then  repeats  the  integration  upward  to  some  z^,  etc. 

If  there  were  no  entrainment  (i.e.,  m  constant),  this  method  would 
have  the  cloud  oscillate  indefinitely.  Actually,  m  increases  with  time 
so  that  the  amplitude  of  oscillation  decreases.  That  is,  entrainment 
exercises  a  braking  effect.  But  the  braking  may  be  insufficient. 
Calculation  with  only  this  form  of  drag  give  excessive  velocities  and 
oscillations.  (Refs.  17,  25,  and  calculations  made  with  the  set  of 
equations  given  in  this  report.) 

The  additional  drag  force  required  can  be  expressed  by  a  virtual 
mass,  a  drag  coefficient,  or  as  a  coefficient  of  momentum  exchange  or 
momentum  loss  by  turbulence.  (Sec.  2.6.2) 
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The  static  and  quasi-static  approaches  give  minimum  and  maximum 
values  of  maximum  cloud  height,  respectively.  Both  are  energy  methods 
and  independent  of  time.  The  first  disregards  momentum  completely. 

The  second  is  an  energy- integral  method,  and  disregards  any  dissipative 
drag  or  friction,  or  either  upward  or  downward  motion. 

B.2  MOMENTUM  EQUATIONS 

When  a  momentum  equation  is  used,  if  the  cloud  is  "heavier  than 
air"  its  rate  of  rise  has  decreased  to  zero,  the  cloud  sinks  until  it 
again  runs  out  of  momentum  (as  it  will  in  a  stable  atmosphere  even 
without  entrainment).  Thus,  the  cloud  performs  damped  oscillations. 

With  sufficient  damping,  no  oscillations  occur  because  zero  rate  of 
rise  is  approached  a symptot i cally . 

If  drag  other  than  that  due  to  entrainment  is  neglected,  this 
approach  gives  the  same  buoyancy-height  history  as  the  quasi-static 
balance  (App.  B.1.2),  but  also  the  time  history  which  quasi-static 
approaches  cannot  provide. 

In  a  momentum  equation,  rate  of  change  of  cloud  momentum  (which 
allows  for  mass  change )  is  balanced  against  buoyancy  force,  and  perhaps 
some  decelerating  force  other  than  mass  change.  Without  the  additional 
drag,  the  calculated  rate  of  cloud  rise  may  be  excessively  high. 

Possibilities  for  the  additional  retardation  term  are: 

1.  Virtual  mass  (potential  flow  theory) 

2.  Drag  coefficient 

3*  Momentum- exchange  or  loss  coefficient. 

B.2.1  Potential  Flow 

The  cloud  is  moving  in  a  fluid  that  is  neither  homogeneous  nor 
infinite,  and  the  cloud  itself  is  fluid,  not  a  rigid  body,  and  there  is 
a  dissipative  loss  of  momentum  to  turbulence  (Sec.  2.6.2).  This  sug¬ 
gests  that  use  of  potential  flow  theory  for  cloud  rise  must  be  modified. 
Nevertheless  the  potential  flow  approach  may  be  useful. 
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In  this  problem  there  is  not  only  accelerated  flow,  but  also 
change  in  virtual  mass, ^  m',  due  to  change  both  in  "sphere"  (cloud) 
volume  and  in  the  ratio  T/Te  *  pe/p,  where  T  and  Te  are,  respectively, 
cloud  and  ambient  temperatures. 

If  these  two  changes  are  considered  then 


d 

dt 


where  u  = 


■|^(m  +  m' ) 

dz 

dt' 


u 


=  V  (pe  -  p)  g 


(B.3) 


Since,  for  flow  around  a  sphere,  virtual  volume  is  one  half  the 
displaced  volume,  i.e.,  m1  =  (m/2)  (pe/p ),  then 


du 

dt 


2 


IT-  Te 

V  +  we 


g  -  2u 


Te  >  1  dm 
T  +  2tf e  m  dt 


Te  \  1  |  T_  dm_  +  m  d(T/Te) 

T  +  2Te  m  |  Te  dt  dt 


(B.4) 


The  last  of  the  three  terms  on  the  right,  which  corresponds  to  the  change 
in  virtual  mass,  was  omitted  in  Eq.  (3l)  of  Ref.  25. 

It,  as  discussed  at  the  end  of  Sec.  2.6.2,  we  allow  only  for  an 
initial  constant  virtual  mass,  m'  •»  iOq'  «*  iDq  T0/2Te0,  then  Eq.  (B.  3) 
becomes 

d1*  m  If  T  _  it  _  u  dm  I  . 

(B.5) 


du 

dt 


m  +  iSq' 


The  effect  of  the  factor  m/(m  +  itiq'  )  in  the  momentum  equation,  Eq.  (3«l), 
of  our  set  of  cloud  rise  equations  (Sec.  3*l)>  is  to  avoid  unreal¬ 
istically  high  initial  acceleration.  In  Eq.  (3*l)  the  presence  of 
water  vapor,  as  well  as  dry  air,  in  the  cloud  and  in  the  atmosphere  is 
taken  into  account,  so  that  the  temperatures  T  and  Ts  in  the  momentum 
equation  are  replaced  by  the  corresponding  virtual  temperatures  T*  and 

m  # 

Te  * 
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The  potential  flow  formula  for  rise  of  a  bubble  with  a  spherical 
leading  face  ^,27  agrees  suprisingly  well  with  observations  of  early 
cloud  rise,  for  yields  up  to  14.5  MT.  At  this  and  higher  yields 
there  may  be  local  transonic  flow  around  equator  of  cloud  and  additional 
pressure  dreg.  (See  Sec.  2.6.4)  The  equation  for  steady- state  rate  of 
rise  of  an  empty  bubble,  of  radius  r,  is 

U  =  2  \  sr  (B.6) 

The  generalization  of  this  to  non-empty  bubbles  is 


If  we  accept  this  rate  of  rise  as  equal  to  the  observed  rate  of  rise, 
then  from  observed  cloud  rise,  one  could  estimate  the  density  ratio 

p/pe- 

For  an  ellipsoidal  cloud,  instead  of  horizontal  radius  r,  the 
radius  of  curvature,  r2/b,  may  be  used,  where  b  is  the  vertical 
radius  of  the  cloud. 

B.2.2  Drag  Coefficients 

Writing  a  momentum  equation  for  the  cloud,  with  a  drag  term,  such 
as  given  by  Eq.  (2.12)  in  Sec  2.6.2, 

-■£-  (mu)  =  V(pe  -  p)  g  -  Cd(|  Pe  u2)rtr2  (B.8) 
where  V  is  the  cloud  volume,  so  that  for  a  sphere 

-ft“(r3pu)  =  r3  (Pe  -  p)  g  -  |  Op  r2  pg  u2  (B-9) 

If  changes  in  r  and  u  are  neglected,  then 

u  =  / 8.  ertpe  ~  Pi  ... 

y  3CD%  (B-10) 
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Note  the  similarity  "between  Eq.  ( 3 . 10 )  and  the  potential-flow  equation, 

dLz 

Eq.  (B.7).  Allowing  for  acceleration  and  maBS  change,  Bince  u  =  — , 

dt 

r(pe  -  p)  8  -  |  C„  pe  u2  =  u2(3p  §  +  r  £|)  +  rpu  g  (B.ll) 

The  two  "maBB  entrainment"  terms  in  parentheses  on  the  right  Bide 
of  Eq.  (B. 11) act  as  additional  drag  terms. 

At  least  for  yields  of  100-KT  and  more,  cloud  diameter  approaches 
or  exceeds  cloud  height,  so  that  in  the  later  stages  of  cloud  rise 
dr _  # 

1/2.  Since  also  then  p  cs  pe  and  Cp  iB  about  0.5, 

31  8  CD  0. 

That  is,  the  drag- coefficient  term  makeB  a  minor  contribution  to  the 
total  drag.  This  is  confirmed  by  Ref.  5*  But  for  high  yields,  early 
in  cloud  rise,  where  we  can  assume  initial  cloud  density  to  be  about 
10$  that  of  ambient  air,  calculated  drag  coefficients  (to  fit  observed 
rate  of  rise)  may  be  greater  than  typical  high-Reynolds-number  values 
for  spheres  of  0.5  or  less.  As  with  the  breakdown  of  potential  flow 
formulas,  one  suspects  the  influence  of  compressibility.  (Sec.  2.6.3) 

In  any  case,  the  UBe  of  the  drag  concept  is  questionable.  The 
cloud  is  a  fluid  body,  and  itB  boundary  layer  is  largely  entrained,  so 
that  although  cloud  momentum  1b  lost  total  kinetic  energy  is  not.  This 
suggests  the  use  of  a  momentum- exchange  or  loss  function  to  describe  the 
kinetic  energy  but  not  the  exact  flow  pattern  of  the  cloud.  Certain 
entrainment  and  eddy-viscosity  ideas  lead  to  the  same  mathematical  form 
of  the  momentum  equation  as  above,  but  with  more  physical- intuitive 
conviction  (Sec.  2.6). 

B.2.3  Momentum  Exchange  and  Momentum  Loss 

B.2.3.1  External  Exchange.  One  possible  momentum- exchange  rule  is  that, 
per  unit  time,  a  certain  fraction  of  cloud  mass  is  exchanged  for  an 
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equal  mass  (or  possibly  an  equal  volume)  of  ambient  air.  Since  the 
ambient  air  is  stationary  there  is  .  net  loss  of  cloud  momentum  per 
mass  (i.e.,  a  dilution  of  momentum)  even  though  mass  is  not  lost. 
Physically:  a  parcel  of  ambient  air  is  entrained  by  the  cloud.  Part  of 
the  parcel  stays  in  the  cloud  to  increase  cloud  mass.  The  rest  of  the 
parcel  is  "mixed  out,"  i.e.  it  returns  to  the  ambient  air,  taking  with 
it  some  of  the  cloud  momentum. 

If  the  "mixed-out"  air  has  been  in  the  cloud  long  enough  to  acquire 
momentum,  should  it  not  acquire  heat  too?  The  answer  is:  not  if  heat 
transfer  requires  more  thorough  mixing.  Picture  the  ambient  parcel  as 
entrained  and  broken  into  two  parts,  and  one  part  as  ejected.  The 
ejected  part  has  gained  momentum,  but  its  interior  is  undisturbed,  so 
no  heat  transfer  has  occurred. 

What  happens  to  the  lost  momentum  represented  by  an  exchange  term 
in  the  momentum  equation?  The  implication  is  that  it  is  eventually 
dissipated  in  the  atmosphere  without  perturbing  cloud  rise. 

The  treatment  of  exchange  factor  as  a  constant  rather  them  velocity- 
dependent  is  appropriate  to  environmental  turbulence.  For  the  rising 
atomic  cloud,  emy  environmental  turbulence  is  negligible  compared  with 
cloud-induced  turbulence  as  em  influence  on  exchange  coefficients.  In 
any  event,  the  exchange  factor  should  be  less  for  large  than  for  small 
clouds  unless  increased  circulation  compensates  for  the  smaller  area- 
to- volume  ratio. 

Q 

A  semi- theoretical  Justification  for  Anderson's0  cloud  rise 
equation  is  as  follows:  in  Sutton's^  diffusion  type  equation  we  take 
the  exponent  1.88  of  altitude,  z,  as  ^  2,  then  we  get 

t  W  /  ,  or  wV4 

a6  in  Ref.  8.  Taking  1.88  »  2  in  Sutton's  equation  means  assuming 
zero  atmospheric  turbulence,  which  is  reasonable  relative  to  cloud 
turbulence . 


69 


B.2.3»2  Internal  Exchange.  Suppose  momentum  exchange  takes  place 
"inside"  the  cloud  rather  than  on  its  surface,  or  that  the  surface 
exchange  1b  with  air  all  of  which  is  promptly  entrained. 

For  surface  exchange,  the  sequence  would  be:  (l)  internal  motion 
causing  momentum  exchange  near  the  cloud  surface,  (2)  circulation  of 
fluid  "lumps"  from  front  and  sides  to  rear  surface  (base)  of  cloud, 

(3)  entrainment  of  more  fluid,  (4)  loss  of  momentum  and  of  kinetic 
energy  of  rise,  and  (5)  reentry  of  the  randomly  moving  lumps  into  the 
cloud. 

Putting  these  ideas  together  suggests  that  the  exchanged  momentum 
and  the  corresponding  kinetic  energy  of  directed  motion  are  not  s i  tiply 
dumped  into  the  atmosphere,  but  go  into  cloud  turbulence.  The  nar¬ 
rowness  of  the  wake  in  nuclear-cloud  photographs  also  suggests  this 
internal  disposal  process. 

The  momentum  is  not  all  transferred  directly  to  turbulent  kinetic 
energy.  Some  of  it  must  go  to  setting  up  the  observed  vortex  circula¬ 
tion.  But  the  vortex  eventually  disappears  with  degradation  of  its 
motion  to  turbulence,  i.e.,  to  smaller,  fluctuating  vortices  or  finally 
to  thermal  energy;  so  in  calculating  the  effect  of  exchanged  or  "lost" 
kinetic  energy  of  rise,  or  accumulated  energy  of  turbulence,  on  the 
horizontal  cloud  diffusion  after  end  of  rise,  we  may  ignore  the  vortex 
formation  and  breakdown  as  intermediate  steps  in  the  energy-conversion 
process. 

The  set  of  cloud- rise  equations  in  Sec.  3  uses  an  internal  exchange 
mechanism.  Since  this  mechanism  is  closely  related  to  the  entrainment 
mechanism,  the  two  processes  are  discussed  together  (Sec.  2.6). 

B.3  ENTRAINMENT 

One  way  to  idealize  entrainment  is  to  consider  it  as  a  snowball 
mechanism,  i.e.  rolling  contact  entrainment.  Then  the  question  arises: 
how  thick  should  the  layer  of  air  that  the  cloud  rolls  up  be?  Pre¬ 
sumably  it  should  be  proportional  to  the  ratio  of  cloud  density  to 
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ambient  density,  as  discussed  in  Sec.  2.6.1.  If,  aside  from  that,  it 
is  constant,  we  get 

1  dV  ^  1 
V  iz  r2 

where  V  cloud  volume,  and  r  •  cloud  radius.  If  thickness  is 
inversely  proportional  to  the  two-way  snowball  curvature,  we  get 

-  —  =  constant  (except  for  the  density  ratio);  i.e., 

V  dz 

Machta's  model. 1 

Another  entrainment  rule  is  that  robe  of  cloud  mass  increase  is 
proportional  to  rate  of  rise  times  cloud  surface  area.**'*  This  rule, 
with  modifications  is  used  in  our  cloud  rise  equations,  in  Sec.  3* 

Since  it  is  closely  related  to  our  model  of  momentum  change,  the  two 
subjects  are  discussed  together  in  Sec.  2.6. 


71 


APPENDIX  C 


PHYSICAL  CONSTANTS  AND  ENVIRONMENTAL  AND  INITIAL 
CONDITIONS  USED  IN  NUMERICAL  SOLUTION  OF  THE  EQUATIONS 

All  quanities  are  in  the  inks  system,  except  explosion  energy  W, 
which  is  in  kilotons  (KT):  1  KT  =  4.18  x  10^  Joules. 

C.l  PHYSICAL  CONSTANTS 

6  -  9.80 

L  =  2.50  x  10^  for  vapor  to  liquid  transition 
L  =  2.83  x  10^  for  vapor  to  solid  transition 
Ra=  287 

6  =  18/29 

C.2  ENVIRONMENTAL  CONDITIONS 

Te  =  Tg,  1_1_  0^  (z  -  for  the  ith  layer  of  a  4- layer 

atmosphere  defined  recursively  by  sea- level  temperature  Teo,  the 
dividing  heights  O,  z^»  Zo,  z^,  and  the  lapse  rates  ct^,  Ofy. 

Tg _ \  e/Ra^i 

Te,  i-l]  f°r  ^he  ith  of  the  4  layers. 

Thus  p  is  defined  recursively.  p0  is  the  only  independent  pressure 
parameter,  besides  those  specifying  Te • 

Teo  and  pQ  are  the  values  at  sea  level,  z  =  0,  and  not  at  initial 
height  zQ,  unless  zQ  =  0.  1 


P  =  Pi.! 
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Note  that  isothermal  layers,  =  0,  can  not  be  used,  but  can  be 
approximated  by  small  absolute  0C,  say  a  =  .0003. 

Xe  =  he  e  es/p,  where  he,  relative  humidity,  is  given  by  a  quadratic 

polynomial  in  z  for 

_  k  h 

z <T  10  and  is  taken  as  0  for  z^>  10 


es> 


C.3 


saturation  vapor  pressure,  is  given  by  the  integrated  Clausius 
Clapeyron  equation: 


e 


s 


6n.  (T/273) 


-5.13 


25.O  (t-273)/t 

e 


TNITIAL  CONDITIONS 


Values  are  assigned  to  V/,  f,  zQ.  These  values,  together  with 
the  environmental  conditions,  determine  the  initial  values  of  m,  x,  V,  r. 
mo  =  “ao  +  "Vo*  where 


y^fW  (4. 18  x  1012) 


mao 


r  T, 


°  -  (T)dT 


~pa' 


T, 


eo 


gfy  (4. 18  x  1012) 

%o  =  rT0 

cpw(T)dT  +  L 


leo 


Cpa  and  Cpy  are  quadratic  polynomials  in  T  whose  coefficients 
are  parameter's  to  be  specified. 

Strictly,  Te(z0)  should  be  used  in  mao  and  m^,  instead  of  Teo>  but 
the  difference  (when  zQ  /  o)  is  unimportant  because  (T0  -  Teo)  is 
large  compared  to  (Te(z0)  -  Teo). 
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The  initial  values  of  x,  V,  and  r  are  given  by 


xo  “ 

VD  =  B®aT0#/p(a o' 
r0  =  OV^/V  )l/3 


Jk 


APPENDIX  D 


GLOSSARY  OF  COMPUTER  PRINTOUT  SYMBOLS 


3ymbol  in  Symbol  in 

printout _ text 


Output  Symbols 


ST  t 
U  u 
X  x 
T  T 
R  r 


Z 

EK 

V 

WT 

TE 

M 

ES 


V 


ra 

es 


Comment 


After  SWITCH  TO  ELLIPSE  R  refers  to 
horizontal  radius,  vertical  radius 
remaining  fixed 


Saturation  water  vapor  pressure  at 
T.  Values  of  ES  printed  for 
T>  373  have  no  physical 
significance 
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Symbol  in 
printout 


Symbol  In 
text 


Comment 


P 

FW 

ED 


partial  pressure  of  water  vapor  In 
cloud 


rate  of  loss  of  kinetic  energy  of 
rise  due  to  eddy  viscosity,  per 
unit  mass. 


Parameter  Symbols 
DST 

K2  kg 

LAMBDA  \ 

Cl,  C2,  C3 


W 

F 

PHI 

TEO 

CHANGE 

DST2 


w 


f 

4 


Leo 


Runge-Kutta  step  size  between 
t  =  0  and  t  =  1. 


Coefficients  of  polynomial  giving 

relative  humidity,  he,  between 

z  =  0,  z  =  10000: 

he  ■  Cl  -  (C2)z  -  (C3)z2 


Sea- level  temperature 

Value  of  ST  after  which  Runge-Kutta 
step  size  changes  to  DST2 
Runge-Kutta  step  size  when 
ST  »  CHANGE 


7 6 


Symbol  in 
printout 


Symbol  in 
text 


Comment 


BO,  Bl,  B2  Coefficients  of  quadratic  polynomial  in 

T  for  cpa. 

DST1  Runge-Kutta  step  size  for  1  ST  CHANGE 

K  Option  in  program.  K  =  1,  3  select  the 

"alternate"  equations,  K  =  1,  2  select 
the  characteristic  length  t  =  r. 

K  =  3>  4  select  Jl  =  100.  This  last 
value  of  £>  in  in  a  dummy  equation  in 
which  some  other  value  of  such  as 
a  Mach- number  correction,  can  be 
inserted. 


Al,  A2,  A3,  A4  ctj_,  0i2>  Q^,  £%  Lapse  rates  in  the  4  layers  of  the 

atmosphere 

LI,  Z2,  Z3  zi,  Z2,  Z3  Dividing  altitudes  of  the  4- layer 

atmosphere 


PO 

ZT 


p0  Sea  level  pressure 

Altitude  such  that  if  Z;>  ZT  (SWITCH 
TO  ELLIPSE)  cloud  vertical  radius 
remains  constant  (nominal  height  of 
tropopause ) 


DO,  Dl,  D2 


Coefficients  of  quadratic  polynomial  in 
T  for  Cpy. 


TF 

C3 


Freezing  point  used  to  select  value  of 
L  (after  SWITCH  TO  WET) 

See  Cl,  C2 
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Symbol  in 
printout 

PRINT 

RK3 


Symbol  in  Comment 

Text 


End  time  of  computation 

If  RK3  =  1  the  momentum  equation 
contains  the  initial  virtual  mass 
factor.  If  RK3  =  0  this  factor  is 
omitted. 
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APPENDIX  E 


COMPUTER  PROGRAM 


The  following  program  (App.  E-2)  has  been  run  on  an  IBM  7094* 

An  IBM  704  version  is  also  available.  The  program  was  written  by 
David  Hutchinson  of  C-E-I-R,  Inc.* 

E.l  COMMENTS  ON  THE  PROGRAM  (By  David  Hutchinson) 

General  Remarks. 

The  integrator  used  is  "Illinois  Library  Routine  D2-DRK1-45-SR" 
"Double  Precision  Floating  Point  Runge  Kutta."  Double  precision  was 
not  needed;  however,  the  routine  was  handy  and  had  been  checked  out  by 
the  programmer  at  the  University  of  Illinois  in  August  1963*  The 
extra  time  required  to  do  the  7094  double-precision  operations  is 
negligible  compared  to  the  time  of  the  DERIV  subroutine  which  is 
single-precision  FORTRAN  coded.  The  equations  are  evaluated  in  double 
precision  (in  the  7094  version). 

In  the  704  version  of  the  integrator  (which  will  run  on  a 
709/90/94),  the  double-precision  7094  instructions  have  been  replaced 
by  the  corresponding  single-precision  instructions.  Now  the  equations 
are  evaluated  in  single-precision.  This  will  cause  a  slight  deteri¬ 
oration  in  accuracy  on  long  runs,  but  not  enough  to  affect  the  number 
of  digits  printed  in  this  problem. 

Because  the  integrator  is  double-precision  and  was  not  designed  to 
be  used  with  FORTRAN,  the  FORTRAN  program  might  look  curious  to  some. 


*Ncw  at  the  University  of  Illinois,  Urbana. 
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All  the  variables  and  their  derivatives  occupy  two  core  locations, 

2n  and  2n  +  1.  The  more  significant  part  is  in  2n  and  the  less 
significant  pert,  which  is  zero  usually,  is  in  2n  +  1.  This  is  so 
because  the  709^  double-precision  instructions  require  the  more 
significant  part  in  an  even  core  location.  Because  FORTRAN  II 
stores  its  arrays  in  decending  core  locations,  the  derivative  for 
example,  has  its  more  significant  part  put  in  DU(2)  by  the  DERIV  sub¬ 
routine.  Hie  less  significant  part  is  zero  which  is  put  into  DU(l) 
once  and  for  all  in  the  MAIN  routine. 
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Punching  the  data  cards. 

All  the  data  are  punche  in  the  column  F  fields  except  K  which  is 
punched  in  a  1  column  I  field  (column  71  on  card  3)*  At  the  Oakland 
installation  of  C.E.I.R.  (where  this  program  was  debugged)  the  version 
II  FORTRAN  monitor  is  used  and  zeros  are  punched  as  0.0.  If  a  blank 
field  is  read  (it  should  be  read  as  minus  zero  for  this  program  to 
run  correctly)  then  the  field  is  ignored  and  the  value  of  the  datum 
which  was  last  read  in  is  used.  Hence,  data  which  do  not  change  from 
one  run  to  the  next  need  not  be  punched  again,  except  for  card  1  where 
all  7  fields  must  be  punched  each  time.  The  five  data  cards  are 
punched  as  follows: 


Cols 

1-10 

11-20 

21-30 

31-40 

41-50 

51-60 

61-70 

71-80 

Card  1 

U 

RK3 

T 

blank 

Z 

EK 

D8T 

Card  P 

K2 

\ 

Cl 

C2 

W 

F 

f* 

Card  3 

TEO 

CHANGE 

DST2 

BO 

B1 

B2 

DST1 

K  (must  be 
in  col  71) 

Card  4 

A1 

A2 

A3 

a4 

Z1 

Z2 

Z3 

PO 

Card  5 

ZT 

DO 

D1 

D2 

TF 

C3 

PRINT 
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DIFFUSION  DECK 


•HUTCHINSON-HUtbSCH  NAVAL  RADIOLOGICAL  LAB. 


« 

FAP 

COUNT 

90 

OLD 

OPD 

044171 100000 

OST 

OPO 

460)71100000 

OF  AD 

OPD 

010171100000 

OF  HP 

OPO 

026171100000 

*  UOI 

DOUBLE 

PRECISION  FL»  PT •  RUNGE-KUT T A-G 1 LL 

ENTRY 

ORK  1 

NOP 

•  REPLACE  THIS  CARO  8Y  EVtN  PSEUDOP  WHEN  IT  BECOMES  OPERATIVE 

002 

DRK1 

SXA 

SR*  1 

SX  A 

SR+1 *2 

NZ1« 

6*4 

TRA 

NEW 

AXT 

8*1 

SXA 

5X4  |4 

EVAL 

^  ALL 

OERIV 

AXT 

**•2 

NLOOP 

OLD 

***2 

OFMP 

## 

OST 

##,2 

13 

OLD 

B+8  *  1 

OFMP 

**•2 

15 

OFAO 

**•2 

016 

OFMP 

A+8  *  1 

OST 

TS 

OFAO 

**•2 

19 

OST 

**  *2 

20 

OLD 

**•2 

021 

OFMP 

C+8.1 

OFAO 

TS 

OFAO 

TS 

OFAO 

TS 

OFAO 

**•2 

26 

OST 

**•2 

27 

T  I  X 

NLOOP *2 *2 

028 

T  I  X 

EVAL* 1*2 

029 

SR 

AXT 

***1 

30 

AXT 

**  *2 

31 

$X4 

AXT 

**•4 

32 

TRA 

7*4 

33 

NEW 

ClA* 

3*4 

34 

STO* 

6*4 

35 

STA 

NLOOP-1 

036 

STA 

ZQ-1 

37 

AOD 

4*4 

38 

STA 

NLOOP+4 

039 

STA 

NLOOP+13 

040 

STA 

NL00P+16 

041 

STA 

ZO 

42 

CLA 

1*4 

43 

ADC* 

5*4 

44 

STA 

NL00P+8 

045 

STA 

NL00P+9 

046 

CLA 

2*4 

47 

ADO* 

5*4 

48 

STA 

NLOOP 

049 

STA 

NLOOP+2 

050 

STA 

NLOOP *3 

031 

STA 

NLOOP *10 

032 
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CLA 

3. A 

53 

STA 

NLOOP+1 

054 

AXT 

**.l 

55 

£0 

ST  2 

**»1 

56 

T1X 

20*1*1 

57 

TRA 

EVAL-2 

#50 

A 

OCT 

200400000000. 145000000000 

59 

OCT 

177453730314*144600610316 

060 

OCT 

201665011714.146637635715 

061 

OCT 

176525252525*143252525253 

062 

B 

OCT 

602400000000.547000000000 

63 

OCT 

60 1400000000. 546000000000 

64 

OCT 

601400000000.546000000000 

65 

OCT 

6024000000UO. 547000000000 

66 

C 

OCT 

600400000000 *545000000000 

67 

OCT 

577453730314.544600610316 

060 

OCT 

601665011714.546637635715 

069 

OCT 

600400000000.545000000000 

70 

TS 

05$ 

2 

71 

END 

72 

* 

FORTRAN 

C 

THIS  IS 

THE  MAIN  ROUTINE 

DIMENSION  PAR  I  37  I *A  I 741*01 16) *DUI 2) *0X12) *DT ( 2 ) *0V ( 2 ) «DZ I  2 ) *DEK( 2  I 
2*SMALLT 12) .DRMI 2 » .DUMMYI 4)  .DWTI2 ) 

COMMON  DUMMY  . A .DWT  .DRM.DU.DX  »DT .DV.D2 .DEK * SMALL T .Q.TE.P.P3.P1.P2.N 

2.TE1.TE2.TE3.ED.NZT.RZT.RK3.0I 
EQUIVALENCE  I DUMMYI 2 ) .WT )  .IDUMMYI4) .RM) 

I  •  IA(2)»U).(A(4).X).IA(6>*T).IA(6)*V).(A(10)«Z).(A(12)*EK 

2» Y ) » I A I IA ) *D5T I  *  I Al 16 ) *RK2 ) • I A  1 18 ) »RL) ♦ ( A  I  20) *C1 ) *  <  A ( 2  2 ) *C2 ) * ( A( 2A 

3)»WI»IA(26)*FI*(A(28)»PHI)» 

Al A I  30) .TEO) • I A I  32  I (CHANGE) *  I A I 3A ) *DST2  ) *  I A I  36 ) . BO ) *  I Al 30 1 *B1» » 

5  IAIAO) *02) »<A(42>»DST1  )•  (AI44) .K) 

6*  I  Al 46 ) *A1 ) *  I A 1 48) *A2 ) *1 Al 50) • A3) • < At  52 ) • A4 ) *  I A 1 54 ) *Zl ) *  I A I  56 ) *22 ) 
7.  I  Al 58 ) *23 ) • (A (60 ) «P0 ) » I AI62 1 *ZT ) « I Al 64 ) .DO) . I A 1 66 ) »D1 ) *  I A 1 66 ) »D2 ) 
8  *  I Al 70 ) *TF) . I A ( 72 ) »C3) *IA|7A) .PRINT  ) 

C  INPUT  THE  INITIAL  CONDITIONS  FOLLOWED  BY  DST  AND  PARAMETERS 

1  READ  INPUT  TAPE  5*19. PAR 

19  FORMAT! 7F 10.0/ 7F 10.0/ 7F 10.0. II /8F10.0/8F10.0) 

C  IF  T*PAR  1 3 )  IS  NEGATIVE  IT  IS  A  SIGNAL  TO  STOP 

I F I  PAR  I  3 )) 18  .4  .4 
18  CALL  DUMP 

A  DO  3  1*1 »37 

C  IF  FIELD  IS  BLANK  USE  LAST  PARAMETER  VALUE  INPUT 

IF  I  PAR  I  I  1)2.33.2 
B  33  D*PARII)+201AOOOCOOOO 
I  FI D)  3.2.2 

2  Al 2*1 ) *PARI I ) 

3  Al 2*1-1 ) *0. 

RK3*X 

C  PRINT  THE  PARAMETERS 

WRITE  OUTPUT  TAPE6.17.IAI I) .1*14.74.2) .RK3 
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17  FORMAT ( 5H1DST“F8.5»9H  K2»F10.7»9H  IAMBDA-F10.6.3H  C1-F11.7.5H 

2C2-F11.7.4H  W“E12.6»4H  F“F12.8»6H  PHI«F12.8/ 

3  7H0  TE0-F10.6#9H  CHANGE“F7.3»7H  DST2-F10. 4»9H  B0-FU.2  » 

A  »H  B1“F11.4»SH  02»E12.4»7H  0ST1-F7.1»*H  K-11/6H0  A1«F8.3» 

33H  A2“F8.5t5H  A3“F8.5»5H  A4“F8.3»3H  Z1“F8.0.3H  Z2«F8.0»3H  Z3 

6“F8.0t5H  PO-F9.0/6HO  ZT-F8.0.9H  D0“F8.1»5H  D1“F8.3»5H  02-F10* 

78  »9H  TF-F8.1.5H  C3-F8.6.8H  PR  I N T “F6« 1 1 6H  RK3-F4. 1/ 1HO) 

B  BN “000000000020 

Fl“0. 

SMALLTI 1 > “0# 

SMALLT ( 2 ) “0* 

DWT-O. 

ORM-O. 

OU-O. 

OX-O. 

DT-O. 

0V“0. 

0Z“0. 

DEK*0. 

TO-T 

DWT (2 ) “0. 

TE“TEO 

TE1“TE0-A1*21 

TE2-TE1-A2MZ2-Z1I 

TE3“TE2-A3*(Z3-Z2> 

XX-9. 8/287. 

P1“PO*(TE1/TEO)**(XX/A1I 
P2“P1*( TE2/TE1 1**1 XX/ A2 I 
P3“P2*(TE3/TE2 I**(kx/A3l 

RMA0“PHI*F*W*4.18E12/CB0*IT0-TE0H>Bl/2.*  I TO*TO-TEO*TEO) 

2  ♦B2/3.  *(TO*TO*TO  -  TEO*TEO*TEOI I 

RMWO“ ( l.-PHI ) *F*W*4. 18E12/I  00  *(TO-TEO»*  01/2. *( TO*TO-TEO*TEO>+ 

2  02/3.  *(T0*T0*T0-TE0*TE0*TE0U2300000.) 

C  SET  INITIAL  CONO.  FOR  X 

X  -  RMWO/RMAO 

C  SET  INITIAL  FOR  V 

P“PO*( (TE0-A1*Z  )/TEO)**(XX/Al I 
V“(RMAO+RMWO )*287»*T0*( 1»+29.*X/18» )/  « P* C l.*X) I 
C  SET  INITIAL  CONO.  FOR  R 

R“ ( 3.* V/ 12. 3663706 1 **.333333333 

RMaRHAO+RHWO 

Q! ■•S*RM*T* ( 18 •♦29»*X I* ( 1»+XE)/(TE*(18»*29. *XE) *( l.+X 1 1 
N-l 
WT"0. 

RZT  — 1.0 

WRITE  OUTPUT  TAPE  6.  99 

99  FORMAT ( 4H  ST»  7X»4HU  »7H  X  »9H  T  »10H  R  » 

2  9H  2  »8H  EX  t  3HV  >8H  WT  »9H  TE  »12H 

3  M  *10H  ES  t  8H  P  tlOH  PW  »9H  £0  /1H  I 

C  PRINT  INITIAL  CONDITIONS 

GO  TO  33 
8  PU-U 

PX“X 
PT“T 
PZ“Z 
PEX“EX 
PV-V 
PRM-RM 
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TAKE  A  RUNGE-KUTTA  STEP 
CALL  DRK1(Y.DEK(2)  »DST.Q(16)  .UN.FLI 

SMALLT  (  2 ) -SMALLT ! 2  I ♦DST 
C  TAKE  SMALL  STEPS  INITIALLY 

IF ( SMALLT (2)-1.0)8»87»88 
•7  DST-DST1 

88  IF(RZT»888.89.89 

89  R«SQRTF(3.#V/ ( RZT*12 . 5663706  I  1 

GO  TO  35 

888  R- I  3. »V/ 12. 5663706  I •». 3333 33333 

35  PW-P*X*29./( 18.  +29.»X) 

ES-611.*  (T/273. I»*(-5.13I»EXPFJ !25.»(T-273. I  I /T I 

11  WRITE  OUTPUT  TAPE  6 *16 .SMALLT 1 2 » .U.X .T .R.Z.EK.  V  »WT  .TE.RM, 
2ES.P.PW.E0 

16  FORMAT  I F6. 1 »F  8. 1  .F9.4 »F8 . 1 »F  8.0.F 8.0.F7 .0 *E 11 • 3 »F8 • 4.F7. 1 .Ell . 3 

2  .F10.0.F9  • 0.F8.O  .F12.3I 

C  N-l.  DRY  MODE  N-2.  WET  MODE  N-3.  SMALL  STEP  DRY  MODE 

IF  (N-2)  150.  ISA.  1531 

150  IF  (ES-PW)  152.152.151 

152  SAVE-DST 
DST-0.5 

C  RESTORE  VARIABLE  VALUES  AT  START  OF  LAST  STEP 

SMALLT I  2 ) "SMALLT 1 2 I-SAVE 
U-PU 
X»PX 
T-PT 
Z-PZ 
EK-PEK 
V-PV 
RM-PRM 
N»3 

C  NOW  TAKE  SMALL  STEPS  UNTIL  ES  LESS  THAN  PW 

GO  TO  8 

1531  IF  I  ES-PW) 41.41*8 
41  DST»SAVE 

N-2 

WRITE  OUTPUT  TAPE  6.77 
77  FORMAT! 14H0SWITCH  TO  WET) 

GO  TO  151 

154  IFIWT+. 00000001)  153.153.151 

153  N-l 
WT-O. 

OWT-O.O 
DWT (2 1-0.0 

WRITE  OUTPUT  TAPE  6*66 
66  FORMAT  I 14H0SWI TCH  TO  ORY) 

C  IF  RZT  IS  POSITIVE  WE  ARE  IN  ELLIPSOIDAL  MODE 

151  IF(RZT) 50.1511. 1511 

C  SWITCH  TO  ELLIPSOIDAL  IF  Z  LARGER  THAN  ZT 

50  IF ( Z~ZT 11511.51.51 

51  RZT-R 

WRITE  OUTPUT  TAPE  6.52 
32  FORMAT! 24H  SWITCH  TO  ELLIPSE.  R-RH) 
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1511  IF<SMALLTI2)-CHANGE)  14  *  15*14 
15  DST  »DST2 

14  I F ( ABSF 111-10*1  1*20*20 

20  IFIR-l.)  1*21*21 

21  IF  (SMALLT 12 )-10.  )22, 210*210 

210  IF(A»SF(DU(2)  )-.l>  211*22*22 

211  I F  t  AOSFI U) - • 1 )  1*22*22 

22  IF  <SMALLTU)-PRINT)13,  13.  1 

13  IFi2-lGOOO.*W«*. 25)8*8*1 

END 

*  fortran 

subroutine  deriv 

C  THIS  PROGRAM  COMPUTES  THE  DERIVATIVES  FOR  DRK1 

C  THIS  PROGRAM  ENTERED  4  TIMES  FOR  EACH  R-K  STEP 

DIMENSION  PAR ( 37 ) *A ( 74 ) *Q! 16 ) »DU I  2 ) »DX! 2) »DT(2 ) *DV(2 ) *0212) *DEK( 2) 
2 .SMALLT ( 2 ) *DRM ( 2 ) .DUMMY ( 4 )  *DWT ( 2 ) 

COMMON  DUMMY  *A *DWT  *DRM*DU*DX*DT  «DV *DZ *DEK .SMALLT *Q»TE  »P»P3»P1 *P2»N 
2.TE1  »TE2*.TE3»ED.NZT»RZT.RK3»OI 

EQUIVALENCE  ( DUMMY ( 2 ) »WT )  *  I  DUMMY  1 4 ) *RM ) 

1  *  (A(2)*U)»(A(4)*X)*(A(6)*T)*(A(8)*V)*(A(10)*Z)*(A(12)*EK 

2* Y  )  * ( A( 14 ) *DST)*(A(16)«RK2)*(A(18)*RL)*(A(20)*C1 ) *  I A ( 22 ) *C2 ) . I A I  24 

3)*W)»(A'26)*F)*(A(28)*PHI)* 

41AI30) *TEO) • f  A ( 32) .CHANGE) »«A( 34) *DST2  ) »(A(36) *80) •  «A( 38) .81) . 

5  I A(  40 ) *82 ) *(A(42) *DST 1  )»  (A(44)*M 

6* ( A ( 46 ) * A 1 ) *(A(48)*A2)*(A(50) * A3 ) ♦ I A( 52 ) *A4) *  (  A ( 54 ) »Z 1 ) *(A(56) *  Z  2 ) 

7»<A(58)»Z3).(A<60) »P0) *( A(62) *  ZT ) * ( A 1641*00) » I  A ( 66 ) *D1 »  * ( A ( 68 i *02) 

8*1 A(70) *TF) *(A(72)*C3)»(A(74) .PRINT) 

DZ 12  )»U 

C  COMPUTE  TE  AND  P 

XX«9. 8/287. 

IF(Z“Z 1)80*80. 81 

80  TE»TE0-A1*Z 

P«PU*( TE/TE0)**(XX/A1 ) 

GO  TO  89 

81  IF(Z~ 22)82*82*83 

82  TE»TE1-A2*IZ-Z1) 

P«D1*(TE/TE1)**(XX/A2) 

GO  <0  89 

83  IF(Z-Z3)84*84»85 

84  TE«TE2-A3*«Z-Z2) 

P»P2*(TE/TE2)**IXX/A3) 

GO  TO  89 

85  TE«TE3-A4*IZ-Z3) 

P»P3*«TE/TE3)**IXX/A4) 

89  CONTINUE 


62  IF(Z-10000. )3C *36.37 

37  XE»0. 

GO  TO  38 

36  XE  »  18.*«C1~C2«Z-C3*Z»Z)* 

1  611*««TE/273* )**I-5.13)*EXPFI I  25 .* ( TE-273. )  )/TE 

2)/(P*29.) 

38  CPAI«BO*I T-TEM-Bl/2.«IT*T-TE*TE)+82/3.*  IT*T*T-TE*TE*TE) 
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CPW*DO  +  D1*T  ♦  D2*T*T 
CP«(bO*bl*T+b2»  T*T  ♦  X*CPW ) / ( X+l • I 

OXE«(  l.+XE)/(l.+29.*XE/18.) 

QX«( l.+29.*X/18. > / 1 1 • ♦X » 

QT-T/TE 

IF ( RZT ) 35*  70.  70 

35  R~( 3. *9/12.5663706) **.333333333 

S9-3./R 
GO  TO  49 

70  R-SQRTF( 3. *V/ ( RZT *  12 .5663706 ) ) 

ECC»SQRTF(R*R-RZT*RZT  >/R  +  1.0E-15 

SV»3.1415926*«2.*R*R+RZT*RZT*  LOGF(  ( l.+ECO/l  l.-ECC)  )/ECC)/V 

49  J«K 

999  07»  MAX1F ( ABSF ( U) .SORTF ( 2»*EK ) )  OIFUSION 

GO  TO  (50 .50.51.51). J 

50  IFIRZT 153.54. 54 

53  R2*R 

GO  TO  52 

54  R2-RZT 
GO  TO  52 

51  R2-100. 

52  GO  TO  (60.61.60.61) »J 

60  DRM(2»«S9*Q7»<1.+29.*X/18.)*T  *RL  *RH*OXE/( ( l.*X*WT)»TE)  DIFUSION 
QQ»1.0 

GO  TO  621 

61  ORM ( 2 ) *  SV*RM*Q7*RL  OIFUSION 

00  «  OT*OX*OXE»( l.+X)/( l.*X*WT) 

621  DU( 2 ) * ( 9.8* (  OT *QX  *  ( 1#*X ) /I l.+X+WT )*OX€  -  1.) 

2  -(00*  07  *2.*RK2/R2  *  DRM(2)/RM)*U  )• 

3  (RM+QI*( 1.-RK3) )/(RM*OI) 

M*N 

GO  TO  ( 100,101.100) .M 

100  OX ( 2 ) --( 1.*X)*IX-XE)*DRM(2)/IRM*( l.+XE) ) 


DT(2)«-0X»0T*9.8»U/CP  *  OXE  -  CPAI *ORH( 2 ) / ( CP*RM ) 

C  RK2  IS  K2 *  RL  IS  LAMBDA.  RM  IS  M 

DV(2)*(9.8*0XE* 

2  ( l./(287.*GX)-l./CP)*QX*U/TE*(l.-CPAI/IT*CP) )*DRM(2)/RM  ♦ 

3  U.*DX(2)/(18.«QX*(1.+X)«(1.+X)  )  )*V 


GO  TO  555 

C  THIS  IS  THE  WET  PART 

101  01  ■  1.  ♦  X*29./18. 

IF(T-TF)102. 103.103 

102  CL*2.83E6 
GO  TO  104 

103  CL*2. 5E6 

104  02  «  CL*X/(  287. *T  ) 

03  ■  18.*r2/(  T*29.  ) 

04  *  1.  *  02 

05  ■  1.+  CL*03/CP 

06  ■  CL*(X-XE!/CP  ♦  T-TE 

DT ( 2 ) *  ( -QX*OT »9.0*O4*U/CP  *QXE  -Q6*DRM < 2 ) /RM ) /G5 
DX ( 2 )  »  01* (03*0T ( 2 )  ♦  9.8*X*U/ ( 287. *TE) *QXE) 

D9(2>»9*(  9.8* ( 1 • / ( 287.*QX ) “  1 ./CP ) *GX*Q**U/ < Q5*TE )  *QXE 
1  +  (1.-  Q6/(T*Q5) )*0RM(2)/RM) 

DWT(2)«-(  1 .+X+WT ) *( WT+X-XE) *DRM( 2 • / ( < I.+XE )*RM)  -  DX<2) 
555  ED*RK2*  U*U*Q7  *00/R2 

DEK ( 2  ) *  ED  -  EX*0RM(2)/RM 
RETURN 
END 
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